Reference documentation for deal.II version Git 73c87d96ef 2021-11-30 22:54:44 +0100
\(\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 
31 namespace 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>
49  : BatchOptimizer()
50  {
51  set_optimization_method(optimization_method, optimization_flags);
52  }
53 
54 
55 
56  template <typename ReturnType>
58  const BatchOptimizer<ReturnType> &other)
59  : method(other.method)
60  , flags(other.flags)
66  {}
67 
68 
69 
70  template <typename ReturnType>
71  void
75  {
76  Assert(
77  optimized() == false,
78  ExcMessage(
79  "Cannot call set_optimization_method() once the optimizer is finalized."));
80 
81 # ifndef HAVE_SYMENGINE_LLVM
82  if (optimization_method == OptimizerType::llvm)
83  {
85  }
86 # endif
89  }
90 
91 
92 
93  template <typename ReturnType>
94  enum OptimizerType
96  {
97  return method;
98  }
99 
100 
101 
102  template <typename ReturnType>
103  enum OptimizationFlags
105  {
106  return flags;
107  }
108 
109 
110 
111  template <typename ReturnType>
112  bool
114  {
116  }
117 
118 
119 
120  template <typename ReturnType>
121  bool
123  {
124  if (dependent_variables_output.size() > 0)
125  {
128  ExcInternalError());
129  return true;
130  }
131 
132  return false;
133  }
134 
135 
136 
137  template <typename ReturnType>
138  bool
140  {
142  }
143 
144 
145 
146  template <typename ReturnType>
147  void
150  {
151  Assert(optimized() == false,
152  ExcMessage(
153  "Cannot register symbols once the optimizer is finalized."));
154 
155 # ifdef DEBUG
156  // Ensure that all of the keys in the map are actually symbolic
157  // in nature
158  for (const auto &entry : substitution_map)
159  {
160  const SD::Expression &symbol = entry.first;
161  Assert(SymEngine::is_a<SymEngine::Symbol>(*(symbol.get_RCP())),
162  ExcMessage("Key entry in map is not a symbol."));
163  }
164 # endif
165  // Merge the two maps, in the process ensuring that there is no
166  // duplication of symbols
167  independent_variables_symbols.insert(substitution_map.begin(),
168  substitution_map.end());
169  }
170 
171 
172 
173  template <typename ReturnType>
174  void
176  const SymEngine::map_basic_basic &substitution_map)
177  {
180  }
181 
182 
183 
184  template <typename ReturnType>
185  void
187  const SD::types::symbol_vector &symbols)
188  {
189  Assert(optimized() == false,
190  ExcMessage(
191  "Cannot register symbols once the optimizer is finalized."));
192 
193  for (const auto &symbol : symbols)
194  {
195  Assert(independent_variables_symbols.find(symbol) ==
197  ExcMessage("Symbol is already in the map."));
199  std::make_pair(symbol, SD::Expression(0.0)));
200  }
201  }
202 
203 
204 
205  template <typename ReturnType>
206  void
208  const SymEngine::vec_basic &symbols)
209  {
212  }
213 
214 
215 
216  template <typename ReturnType>
219  {
221  }
222 
223 
224 
225  template <typename ReturnType>
226  std::size_t
228  {
229  return independent_variables_symbols.size();
230  }
231 
232 
233 
234  template <typename ReturnType>
235  void
237  {
238  Assert(optimized() == false,
239  ExcMessage(
240  "Cannot register functions once the optimizer is finalized."));
241 
242  register_scalar_function(function);
243  }
244 
245 
246 
247  template <typename ReturnType>
248  void
251  {
252  Assert(optimized() == false,
253  ExcMessage(
254  "Cannot register functions once the optimizer is finalized."));
255 
256  register_vector_functions(functions);
257  }
258 
259 
260 
261  template <typename ReturnType>
262  void
264  const SymEngine::vec_basic &functions)
265  {
268  }
269 
270 
271 
272  template <typename ReturnType>
275  {
277  }
278 
279 
280 
281  template <typename ReturnType>
282  std::size_t
284  {
285  if (has_been_serialized == false)
286  {
287  // If we've had to augment our map after serialization, then
288  // this check, unfortunately, cannot be performed.
291  ExcInternalError());
292  }
293  return dependent_variables_functions.size();
294  }
295 
296 
297 
298  template <typename ReturnType>
299  void
301  {
302  Assert(optimized() == false,
303  ExcMessage("Cannot call optimize() more than once."));
304 
305  // Create and configure the optimizer
308 
309  const SD::types::symbol_vector symbol_vec =
312  *opt = dynamic_cast<typename internal::DictionaryOptimizer<
313  ReturnType>::OptimizerType *>(optimizer.get()))
314  {
316  ExcInternalError());
317  internal::OptimizerHelper<ReturnType,
319  initialize(opt,
321  symbol_vec),
325  }
327  *opt = dynamic_cast<typename internal::LambdaOptimizer<
328  ReturnType>::OptimizerType *>(optimizer.get()))
329  {
331  ExcInternalError());
332  internal::OptimizerHelper<ReturnType,
334  initialize(opt,
336  symbol_vec),
340  }
341 # ifdef HAVE_SYMENGINE_LLVM
343  *opt = dynamic_cast<typename internal::LLVMOptimizer<
344  ReturnType>::OptimizerType *>(optimizer.get()))
345  {
347  ExcInternalError());
348  internal::OptimizerHelper<ReturnType,
349  internal::LLVMOptimizer<ReturnType>>::
350  initialize(opt,
352  symbol_vec),
356  }
357 # endif
358  else
359  {
360  AssertThrow(false, ExcMessage("Unknown optimizer type."));
361  }
362 
363  // The size of the outputs is now fixed, as is the number and
364  // order of the symbols to be substituted.
365  // Note: When no optimisation is actually used (i.e. optimization_method()
366  // == off and use_symbolic_CSE() == false), we could conceptually go
367  // without this data structure. However, since the user expects to perform
368  // substitution of all dependent variables in one go, we still require it
369  // for intermediate storage of results.
371  }
372 
373 
374 
375  template <typename ReturnType>
376  void
379  {
380  Assert(
381  optimized() == true,
382  ExcMessage(
383  "The optimizer is not configured to perform substitution. "
384  "This action can only performed after optimize() has been called."));
386 
387  // Check that the registered symbol map and the input map are compatible
388  // with one another
389 # ifdef DEBUG
390  const SD::types::symbol_vector symbol_sub_vec =
391  Utilities::extract_symbols(substitution_map);
392  const SD::types::symbol_vector symbol_vec =
394  Assert(symbol_sub_vec.size() == symbol_vec.size(),
395  ExcDimensionMismatch(symbol_sub_vec.size(), symbol_vec.size()));
396  for (unsigned int i = 0; i < symbol_sub_vec.size(); ++i)
397  {
398  Assert(numbers::values_are_equal(symbol_sub_vec[i], symbol_vec[i]),
399  ExcMessage(
400  "The input substitution map is either incomplete, or does "
401  "not match that used in the register_symbols() call."));
402  }
403 # endif
404 
405  // Extract the values from the substitution map, and use the other
406  // function
407  const std::vector<ReturnType> values =
408  Utilities::extract_values<ReturnType>(substitution_map);
409  substitute(values);
410  }
411 
412 
413 
414  template <typename ReturnType>
415  void
417  const SymEngine::map_basic_basic &substitution_map) const
418  {
419  substitute(
421  }
422 
423 
424 
425  template <typename ReturnType>
426  void
428  const SD::types::symbol_vector &symbols,
429  const std::vector<ReturnType> & values) const
430  {
431  // Zip the two vectors and use the other function call
432  // This ensures the ordering of the input vectors matches that of the
433  // stored map.
434  substitute(make_substitution_map(symbols, values));
435  }
436 
437 
438 
439  template <typename ReturnType>
440  void
442  const SymEngine::vec_basic & symbols,
443  const std::vector<ReturnType> &values) const
444  {
446  symbols),
447  values);
448  }
449 
450 
451 
452  template <typename ReturnType>
453  void
455  const std::vector<ReturnType> &substitution_values) const
456  {
457  Assert(
458  optimized() == true,
459  ExcMessage(
460  "The optimizer is not configured to perform substitution. "
461  "This action can only performed after optimize() has been called."));
463  Assert(substitution_values.size() == independent_variables_symbols.size(),
464  ExcDimensionMismatch(substitution_values.size(),
466 
468  *opt = dynamic_cast<typename internal::DictionaryOptimizer<
469  ReturnType>::OptimizerType *>(optimizer.get()))
470  {
472  ExcInternalError());
473  internal::OptimizerHelper<ReturnType,
475  substitute(opt, dependent_variables_output, substitution_values);
476  }
478  *opt = dynamic_cast<typename internal::LambdaOptimizer<
479  ReturnType>::OptimizerType *>(optimizer.get()))
480  {
482  ExcInternalError());
483  internal::OptimizerHelper<ReturnType,
485  substitute(opt, dependent_variables_output, substitution_values);
486  }
487 # ifdef HAVE_SYMENGINE_LLVM
489  *opt = dynamic_cast<typename internal::LLVMOptimizer<
490  ReturnType>::OptimizerType *>(optimizer.get()))
491  {
493  ExcInternalError());
494  internal::OptimizerHelper<ReturnType,
495  internal::LLVMOptimizer<ReturnType>>::
496  substitute(opt, dependent_variables_output, substitution_values);
497  }
498 # endif
499  else
500  {
501  AssertThrow(false, ExcNotImplemented());
502  }
503 
505  }
506 
507 
508 
509  template <typename ReturnType>
510  const std::vector<ReturnType> &
512  {
513  Assert(
514  values_substituted() == true,
515  ExcMessage(
516  "The optimizer is not configured to perform evaluation. "
517  "This action can only performed after substitute() has been called."));
518 
520  }
521 
522 
523 
524  template <typename ReturnType>
525  ReturnType
527  const Expression & func,
528  const std::vector<ReturnType> &cached_evaluation) const
529  {
530  Assert(
531  values_substituted() == true,
532  ExcMessage(
533  "The optimizer is not configured to perform evaluation. "
534  "This action can only performed after substitute() has been called."));
535 
536  // TODO[JPP]: Find a way to fix this bug that crops up in serialization
537  // cases, e.g. symengine/batch_optimizer_05. Even though the entry is
538  // in the map, it can only be found by an exhaustive search and string
539  // comparison. Why? Because the leading zero coefficient may seemingly
540  // be dropped (or added) at any time.
541  //
542  // Just this should theoretically work:
543  const typename map_dependent_expression_to_vector_entry_t::const_iterator
544  it = map_dep_expr_vec_entry.find(func);
545 
546  // But instead we are forced to live with this abomination, and its
547  // knock-on effects:
548  if (has_been_serialized && it == map_dep_expr_vec_entry.end())
549  {
550  // Some SymEngine operations might return results with a zero leading
551  // coefficient. Upon serialization, this might be dropped, meaning
552  // that when we reload the expressions they now look somewhat
553  // different to as before. If all data that the user uses is
554  // guaranteed to either have been serialized or never serialized, then
555  // there would be no problem. However, users might rebuild their
556  // dependent expression and just reload the optimizer. This is
557  // completely legitimate. But in this scenario we might be out of sync
558  // with the expressions. This is not great. So we take the nuclear
559  // approach, and run everything through a serialization operation to
560  // see if we can homogenize all of the expressions such that they look
561  // the same in string form.
562  auto serialize_and_deserialize_expression =
563  [](const Expression &old_expr) {
564  std::ostringstream oss;
565  {
566  boost::archive::text_oarchive oa(oss,
567  boost::archive::no_header);
568  oa << old_expr;
569  }
570 
571  Expression new_expr;
572  {
573  std::istringstream iss(oss.str());
574  boost::archive::text_iarchive ia(iss,
575  boost::archive::no_header);
576 
577  ia >> new_expr;
578  }
579 
580  return new_expr;
581  };
582 
583  const Expression new_func =
584  serialize_and_deserialize_expression(func);
585 
586  // Find this in the map, while also making sure to compactify all map
587  // entries. If we find the entry that we're looking for, then we
588  // (re-)add the input expression into the map, and do the proper
589  // search again. We should only need to do this once per invalid
590  // entry, as the corrected entry is then cached in the map.
591  for (const auto &e : map_dep_expr_vec_entry)
592  {
593  const Expression new_map_expr =
594  serialize_and_deserialize_expression(e.first);
595 
596  // Add a new map entry and re-search. This is guaranteed to
597  // return a valid entry. Note that we must do a string comparison,
598  // because the data structures that form the expressions might
599  // still be different.
600  if (new_func.get_value().__str__() ==
601  new_map_expr.get_value().__str__())
602  {
603  map_dep_expr_vec_entry[func] = e.second;
604  return evaluate(func);
605  }
606  }
607 
608  AssertThrow(
609  false,
610  ExcMessage(
611  "Still cannot find map entry, and there's no hope to recover from this situation."));
612  }
613 
614  Assert(it != map_dep_expr_vec_entry.end(),
615  ExcMessage("Function has not been registered."));
616  Assert(it->second < n_dependent_variables(), ExcInternalError());
617 
618  return cached_evaluation[it->second];
619  }
620 
621 
622 
623  template <typename ReturnType>
624  ReturnType
626  {
627  return extract(func, dependent_variables_output);
628  }
629 
630 
631 
632  template <typename ReturnType>
633  std::vector<ReturnType>
635  const std::vector<Expression> &funcs,
636  const std::vector<ReturnType> &cached_evaluation) const
637  {
638  std::vector<ReturnType> out;
639  out.reserve(funcs.size());
640 
641  for (const auto &func : funcs)
642  out.emplace_back(extract(func, cached_evaluation));
643 
644  return out;
645  }
646 
647 
648 
649  template <typename ReturnType>
650  std::vector<ReturnType>
652  const std::vector<Expression> &funcs) const
653  {
654  return extract(funcs, dependent_variables_output);
655  }
656 
657 
658 
659  template <typename ReturnType>
660  bool
662  const SD::Expression &func) const
663  {
665  }
666 
667 
668 
669  template <typename ReturnType>
670  bool
672  const SymEngine::RCP<const SymEngine::Basic> &func) const
673  {
674  // SymEngine's internal constants are the valid
675  // reusable return types for various derivative operations
676  // See
677  // https://github.com/symengine/symengine/blob/master/symengine/constants.h
678  if (SymEngine::is_a<SymEngine::Constant>(*func))
679  return true;
680  if (&*func == &*SymEngine::zero)
681  return true;
682  if (&*func == &*SymEngine::one)
683  return true;
684  if (&*func == &*SymEngine::minus_one)
685  return true;
686  if (&*func == &*SymEngine::I)
687  return true;
688  if (&*func == &*SymEngine::Inf)
689  return true;
690  if (&*func == &*SymEngine::NegInf)
691  return true;
692  if (&*func == &*SymEngine::ComplexInf)
693  return true;
694  if (&*func == &*SymEngine::Nan)
695  return true;
696 
697  return false;
698  }
699 
700 
701 
702  template <typename ReturnType>
703  void
705  const SD::Expression &func)
706  {
707  Assert(
708  dependent_variables_output.size() == 0,
709  ExcMessage(
710  "Cannot register function as the optimizer has already been finalized."));
712  const bool entry_registered =
713  (map_dep_expr_vec_entry.find(func) != map_dep_expr_vec_entry.end());
714 # ifdef DEBUG
715  if (entry_registered == true &&
717  Assert(entry_registered,
718  ExcMessage("Function has already been registered."));
719 # endif
720  if (entry_registered == false)
721  {
722  dependent_variables_functions.push_back(func);
723  map_dep_expr_vec_entry[func] =
725  }
726  }
727 
728 
729 
730  template <typename ReturnType>
731  void
733  const SD::types::symbol_vector &funcs)
734  {
735  Assert(
736  dependent_variables_output.size() == 0,
737  ExcMessage(
738  "Cannot register function as the optimizer has already been finalized."));
739  const std::size_t n_dependents_old = n_dependent_variables();
740  dependent_variables_output.reserve(n_dependents_old + funcs.size());
741  dependent_variables_functions.reserve(n_dependents_old + funcs.size());
742 
743  for (const auto &func : funcs)
744  {
745  const bool entry_registered =
746  (map_dep_expr_vec_entry.find(func) != map_dep_expr_vec_entry.end());
747 # ifdef DEBUG
748  if (entry_registered == true &&
750  Assert(entry_registered,
751  ExcMessage("Function has already been registered."));
752 # endif
753  if (entry_registered == false)
754  {
755  dependent_variables_functions.push_back(func);
756  map_dep_expr_vec_entry[func] =
758  }
759  }
760  }
761 
762 
763 
764  template <typename ReturnType>
765  void
767  std::unique_ptr<SymEngine::Visitor> &optimizer)
768  {
769  Assert(!optimizer, ExcMessage("Optimizer has already been created."));
770 
773  {
774  using Optimizer_t =
776  optimizer.reset(new Optimizer_t());
777  }
779  {
780  using Optimizer_t =
782  optimizer.reset(new Optimizer_t());
783  }
785  {
786 # ifdef HAVE_SYMENGINE_LLVM
787  if (internal::LLVMOptimizer<ReturnType>::supported_by_LLVM)
788  {
789  using Optimizer_t =
791  optimizer.reset(new Optimizer_t());
792  }
793  else
794  {
796  }
797 # else
799 # endif
800  }
801  else
802  {
803  AssertThrow(false, ExcMessage("Unknown optimizer selected."));
804  }
805  }
806 
807  } // namespace SD
808 } // namespace Differentiation
809 
810 
811 /* --- Explicit instantiations --- */
812 # include "symengine_optimizer.inst"
813 
814 
816 
817 #endif // DEAL_II_WITH_SYMENGINE
const std::vector< ReturnType > & evaluate() const
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:925
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
bool use_symbolic_CSE(const enum OptimizationFlags &flags)
SD::types::symbol_vector convert_basic_vector_to_expression_vector(const SymEngine::vec_basic &symbol_vector)
const types::symbol_vector & get_dependent_functions() const
void register_scalar_function(const SD::Expression &function)
static const types::blas_int one
static ::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
enum OptimizationFlags optimization_flags() const
types::symbol_vector get_independent_symbols() const
const SymEngine::RCP< const SymEngine::Basic > & get_RCP() const
void register_functions(const types::symbol_vector &functions)
std::vector< ReturnType > dependent_variables_output
static ::ExceptionBase & ExcMessage(std::string arg1)
enum OptimizerType optimization_method() const
#define Assert(cond, exc)
Definition: exceptions.h:1461
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void register_function(const Expression &function)
void substitute(const types::substitution_map &substitution_map) const
static ::ExceptionBase & ExcSymEngineLLVMReturnTypeNotSupported()
const SymEngine::Basic & get_value() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
SymEngine::vec_basic convert_expression_vector_to_basic_vector(const SD::types::symbol_vector &symbol_vector)
map_dependent_expression_to_vector_entry_t map_dep_expr_vec_entry
std::vector< SD::Expression > symbol_vector
SD::types::symbol_vector extract_symbols(const SD::types::substitution_map &substitution_values)
types::substitution_map independent_variables_symbols
ReturnType extract(const Expression &func, const std::vector< ReturnType > &cached_evaluation) const
types::substitution_map make_substitution_map(const Expression &symbol, const Expression &value)
static ::ExceptionBase & ExcSymEngineLLVMNotAvailable()
void register_vector_functions(const types::symbol_vector &functions)
void create_optimizer(std::unique_ptr< SymEngine::Visitor > &optimizer)
std::unique_ptr< SymEngine::Visitor > optimizer
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
bool is_valid_nonunique_dependent_variable(const SD::Expression &function) const
SD::types::substitution_map convert_basic_map_to_expression_map(const SymEngine::map_basic_basic &substitution_map)
void register_symbols(const types::substitution_map &substitution_map)
void set_optimization_method(const enum OptimizerType &optimization_method, const enum OptimizationFlags &optimization_flags=OptimizationFlags::optimize_all)
static ::ExceptionBase & ExcNotImplemented()
static const types::blas_int zero
std::map< SD::Expression, SD::Expression, internal::ExpressionKeyLess > substitution_map
types::symbol_vector dependent_variables_functions
int(&) functions(const void *v1, const void *v2)
static ::ExceptionBase & ExcInternalError()