Small tweaks

doc/tut_expression_system.rst

@@ -19,14 +19,6 @@ Constants and parameters are mathematically equivalent, the only difference bein
 that the value of a constant is determined when the expression is created, whereas
 the value of a parameter is loaded from a user-supplied data array at a later stage.
 
-As a simple example, here is a graphical representation of the
-AST for the expression :math:`\left( 1 - x^2\right)y-x`:
-
-.. image:: images/ast.png
-  :scale: 30%
-  :align: center
-  :alt: Example of a symbolic tree in heyoka
-
 Thanks to the use of modern C++ features (user-defined literals, overloaded operators, etc.),
 the construction of the AST of an expression in heyoka can be accomplished via natural
 mathematical notation:

include/heyoka/expression.hpp

@@ -526,7 +526,7 @@ HEYOKA_BEGIN_NAMESPACE
 
 // NOTE: when documenting, we need to point out that the expressions
 // returned by this function are optimised for evaluation. The users
-// can always normalise() these expressions if needed.
+// can always unfix() and normalise() these expressions if needed.
 template <typename... KwArgs>
 dtens diff_tensors(const std::vector<expression> &v_ex, KwArgs &&...kw_args)
 {

src/expression_diff.cpp

@@ -1073,43 +1073,22 @@ auto diff_tensors_impl(const std::vector<expression> &v_ex, const std::vector<ex
 
     get_logger()->trace("dtens creation runtime: {}", sw);
 
-    sw.reset();
-
-    // Unfix the derivatives.
-    //
-    // NOTE: unfixing here is mostly a matter of aestethics. At the moment,
-    // I don't see any reason why we would want to keep the expressions fixed
-    // on output, as 1) if these derivatives are iterated to higher orders,
+    // NOTE: it is unclear at this time if it makes sense here to unfix() the derivatives.
+    // This would have only a cosmetic purpose, as 1) if these derivatives are iterated to higher orders,
     // during decomposition we will be unfixing anyway and 2) unfixing does not
     // imply normalisation, thus we are not changing the symbolic structure
     // of the expressions and pessimising their evaluation.
     //
-    // NOTE: keep also in mind that at this time the derivatives are not *fully* fixed,
-    // in the sense that the entries in subs_map are not themselves fixed, and thus when we do
-    // the substition to create the final expressions, we have inserting unfixed
-    // expressions in the result. This ultimately does not matter as subs() does not
-    // do any normalisation by default, but it would matter if we wanted to manipulate
-    // further the expression of the derivatives while relying on fix() being applied
-    // correctly everywhere.
-    std::vector<expression> unfix_diffs;
-    unfix_diffs.reserve(diff_map.size());
-    for (auto &p : diff_map) {
-        unfix_diffs.push_back(std::move(p.second));
-    }
-    unfix_diffs = unfix(unfix_diffs);
-
-    // Construct the return value.
-    dtens_map_t::sequence_type retval_seq;
-    retval_seq.reserve(unfix_diffs.size());
-    auto diff_map_it = diff_map.begin();
-    for (decltype(unfix_diffs.size()) i = 0; i < unfix_diffs.size(); ++i, ++diff_map_it) {
-        retval_seq.emplace_back(diff_map_it->first, std::move(unfix_diffs[i]));
-    }
-    assert(diff_map_it == diff_map.end());
-    dtens_map_t retval;
-    retval.adopt_sequence(boost::container::ordered_unique_range_t{}, std::move(retval_seq));
-
-    get_logger()->trace("dtens unfix runtime: {}", sw);
+    // Keep also in mind that at this time the derivatives are not
+    // *fully* fixed, in the sense that the entries in subs_map are not themselves fixed, and thus when we do the
+    // substition to create the final expressions, we have inserting unfixed expressions in the result. This ultimately
+    // does not matter as subs() does not do any normalisation by default, but it would matter if we wanted to
+    // manipulate further the expression of the derivatives while relying on fix() being applied correctly everywhere.
+    // In any case, if we decide to unfix() here at some point we should also consider unfixing
+    // the expressions returned by models, for consistency.
+
+    // Assemble and return the result.
+    dtens_map_t retval(boost::container::ordered_unique_range_t{}, diff_map.begin(), diff_map.end());
 
     // Check sorting.
     assert(std::is_sorted(retval.begin(), retval.end(),

src/math/prod.cpp

@@ -1016,7 +1016,7 @@ expression prod_to_div_impl(funcptr_map<expression> &func_map, const expression
                     return it->second;
                 }
 
-                // Recursively transform product into divisions
+                // Recursively transform products into divisions
                 // in the arguments.
                 std::vector<expression> new_args;
                 new_args.reserve(v.args().size());

test/model_nbody.cpp

@@ -7,12 +7,16 @@
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #include <stdexcept>
+#include <string>
 #include <variant>
 #include <vector>
 
+#include <boost/algorithm/string.hpp>
+
 #include <xtensor/xadapt.hpp>
 #include <xtensor/xview.hpp>
 
+#include <heyoka/detail/sub.hpp>
 #include <heyoka/detail/sum_sq.hpp>
 #include <heyoka/expression.hpp>
 #include <heyoka/func.hpp>
@@ -60,23 +64,32 @@ TEST_CASE("nbody")
 
         auto ta = heyoka::taylor_adaptive{dyn, n_ic, kw::compact_mode = true};
 
+        // Check that llvm.pow appears only 3 times: its declaration plus 2 uses
+        // for determining the timestep size.
+        std::vector<boost::iterator_range<std::string::const_iterator>> pow_matches;
+        boost::find_all(pow_matches, ta.get_llvm_state().get_ir(), "@llvm.pow");
+        REQUIRE(pow_matches.size() == 3u);
+
         llvm_state s;
         std::vector<expression> vars;
         for (const auto &p : dyn) {
             vars.push_back(p.first);
         }
 
-        // Check that all sums were replaced by sums of squares.
-        auto n_sums = 0, n_sum_sqs = 0;
+        // Check that all sums were replaced by sums of squares, and check the
+        // number of subtractions.
+        auto n_sums = 0, n_sum_sqs = 0, n_subs = 0;
         for (const auto &[s_ex, _] : ta.get_decomposition()) {
             if (const auto *fptr = std::get_if<func>(&s_ex.value())) {
                 n_sums += static_cast<int>(fptr->extract<detail::sum_impl>() != nullptr);
                 n_sum_sqs += static_cast<int>(fptr->extract<detail::sum_sq_impl>() != nullptr);
+                n_subs += static_cast<int>(fptr->extract<detail::sub_impl>() != nullptr);
             }
         }
 
         REQUIRE(n_sum_sqs == 15);
         REQUIRE(n_sums == 18);
+        REQUIRE(n_subs == 45);
         REQUIRE(ta.get_decomposition().size() == 270u);
 
         const auto dc = add_cfunc<double>(s, "cf", {en_ex}, kw::vars = vars);