Added caching for substitute rules. Made substitute args pass-by-refe…

…rence.

core/algorithms/substitute.cc

@@ -13,6 +13,11 @@
 
 using namespace cadabra;
 
+
+/* Global instance of substitute::rules */
+substitute::Rules replacement_rules;
+
+
 substitute::substitute(const Kernel& k, Ex& tr, Ex& args_, bool partial)
 	: Algorithm(k, tr), comparator(k.properties), args(args_), sort_product_(k, tr), partial(partial)
 	{
@@ -21,68 +26,81 @@ substitute::substitute(const Kernel& k, Ex& tr, Ex& args_, bool partial)
 
 	// Stopwatch sw;
 	// sw.start();
-	cadabra::do_list(args, args.begin(), [&](Ex::iterator arrow) {
-		//args.print_recursive_treeform(std::cerr, arrow);
-		if(*arrow->name!="\\arrow" && *arrow->name!="\\equals")
-			throw ArgumentException("substitute: Argument is neither a replacement rule nor an equality.");
-
-		sibling_iterator lhs=args.begin(arrow);
-		sibling_iterator rhs=lhs;
-		rhs.skip_children();
-		++rhs;
-
-		if(*lhs->name=="") { // replacing a sub or superscript
-			lhs=tr.flatten_and_erase(lhs);
-			}
-		if(*rhs->name=="") { // replacing with a sub or superscript
-			rhs=tr.flatten_and_erase(rhs);
-			}
 
-		try {
-			if(*lhs->multiplier!=1) {
-				throw ArgumentException("substitute: No numerical pre-factors allowed on lhs of replacement rule.");
+	// Check if args are present in global_rules
+	// bool skipchecks = k.replacement_rules->is_present(args);
+	bool skipchecks = replacement_rules.is_present(args);
+
+	// skip if args have already been checked
+	if (!skipchecks) {
+		cadabra::do_list(args, args.begin(), [&](Ex::iterator arrow) {
+			//args.print_recursive_treeform(std::cerr, arrow);
+			if(*arrow->name!="\\arrow" && *arrow->name!="\\equals")
+				throw ArgumentException("substitute: Argument is neither a replacement rule nor an equality.");
+
+			sibling_iterator lhs=args.begin(arrow);
+			sibling_iterator rhs=lhs;
+			rhs.skip_children();
+			++rhs;
+
+			if(*lhs->name=="") { // replacing a sub or superscript
+				lhs=tr.flatten_and_erase(lhs);
+				}
+			if(*rhs->name=="") { // replacing with a sub or superscript
+				rhs=tr.flatten_and_erase(rhs);
 				}
-			// test validity of lhs and rhs
-			iterator lhsit=lhs, stopit=lhs;
-			stopit.skip_children();
-			++stopit;
-			while(lhsit!=stopit) {
-				if(lhsit->is_object_wildcard()) {
-					if(tr.number_of_children(lhsit)>0) {
-						throw ArgumentException("substitute: Object wildcards cannot have child nodes.");
+
+			try {
+				if(*lhs->multiplier!=1) {
+					throw ArgumentException("substitute: No numerical pre-factors allowed on lhs of replacement rule.");
+					}
+				// test validity of lhs and rhs
+				iterator lhsit=lhs, stopit=lhs;
+				stopit.skip_children();
+				++stopit;
+				while(lhsit!=stopit) {
+					if(lhsit->is_object_wildcard()) {
+						if(tr.number_of_children(lhsit)>0) {
+							throw ArgumentException("substitute: Object wildcards cannot have child nodes.");
+							}
 						}
+					++lhsit;
 					}
-				++lhsit;
-				}
-			lhsit=rhs;
-			stopit=rhs;
-			stopit.skip_children();
-			++stopit;
-			while(lhsit!=stopit) {
-				if(lhsit->is_object_wildcard()) {
-					if(tr.number_of_children(lhsit)>0) {
-						throw ArgumentException("substitute: Object wildcards cannot have child nodes.");
+				lhsit=rhs;
+				stopit=rhs;
+				stopit.skip_children();
+				++stopit;
+				while(lhsit!=stopit) {
+					if(lhsit->is_object_wildcard()) {
+						if(tr.number_of_children(lhsit)>0) {
+							throw ArgumentException("substitute: Object wildcards cannot have child nodes.");
+							}
 						}
+					++lhsit;
 					}
-				++lhsit;
-				}
 
-			// check whether there are dummies.
-			index_map_t ind_free, ind_dummy;
-			classify_indices(lhs, ind_free, ind_dummy);
-			lhs_contains_dummies[arrow]= ind_dummy.size()>0;
-			ind_free.clear();
-			ind_dummy.clear();
-			if(rhs!=tr.end()) {
-				classify_indices(rhs, ind_free, ind_dummy);
-				rhs_contains_dummies[arrow]=ind_dummy.size()>0;
+				// check whether there are dummies.
+				index_map_t ind_free, ind_dummy;
+				classify_indices(lhs, ind_free, ind_dummy);
+				lhs_contains_dummies[arrow]= ind_dummy.size()>0;
+				ind_free.clear();
+				ind_dummy.clear();
+				if(rhs!=tr.end()) {
+					classify_indices(rhs, ind_free, ind_dummy);
+					rhs_contains_dummies[arrow]=ind_dummy.size()>0;
+					}
 				}
-			}
-		catch(std::exception& er) {
-			throw ArgumentException(std::string("substitute: Index error in replacement rule. ")+er.what());
-			}
-		return true;
-		});
+			catch(std::exception& er) {
+				throw ArgumentException(std::string("substitute: Index error in replacement rule. ")+er.what());
+				}
+			return true;
+			});
+		replacement_rules.store(args, lhs_contains_dummies, rhs_contains_dummies);
+		}
+	else {
+		replacement_rules.retrieve(args, lhs_contains_dummies, rhs_contains_dummies);
+		}
+
 	// sw.stop();
 	// std::cerr << "preparation took " << sw << std::endl;
 	}
@@ -436,3 +454,85 @@ Algorithm::result_t substitute::apply(iterator& st)
 	return result_t::l_applied;
 	}
 
+
+
+
+
+void substitute::Rules::store(Ex& rules,
+								std::map<iterator, bool>& lhs_contains_dummies,
+								std::map<iterator, bool>& rhs_contains_dummies) {
+
+	// if number of stored rules has grown large, clean them up.
+	if (properties.size() >= cleanup_threshold) {
+		cleanup();
+	}
+	// If that didn't fix it, double the cleanup_threshold up to max_size
+	if (cleanup_threshold != max_size && properties.size() >= cleanup_threshold) {
+			if (cleanup_threshold * 2 < max_size) {
+				cleanup_threshold *= 2;
+				} 
+			else {
+				cleanup_threshold = max_size;
+			}
+		}
+	// If we're too big, don't add anything else.
+	if (properties.size() >= max_size) {
+		return;
+	}
+
+	std::weak_ptr<Ex> rules_ptr = rules.shared_from_this();
+	properties[rules_ptr] = { lhs_contains_dummies, rhs_contains_dummies };
+	// Set state of rules to l_checkpointed to track if the rules ever change
+	rules.reset_state();
+	}
+
+void substitute::Rules::retrieve(Ex& rules,
+								std::map<iterator, bool>& lhs_contains_dummies,
+								std::map<iterator, bool>& rhs_contains_dummies) {
+
+	// Rules::present is assumed to have been called to check that the rules are valid
+	std::weak_ptr<Ex> rules_ptr = rules.shared_from_this();
+	lhs_contains_dummies = properties[rules_ptr].first;
+	rhs_contains_dummies = properties[rules_ptr].second;
+	}
+
+
+bool substitute::Rules::is_present(Ex& rules) {
+	// Look to see if the rules are present in the map
+
+	std::weak_ptr<Ex> rules_ptr = rules.shared_from_this();
+	bool rule_found = (properties.find(rules_ptr) != properties.end());
+	if (!rule_found) return false;
+
+	// rules should have l_checkpointed set
+	bool rule_unchanged = (rules.state() == result_t::l_checkpointed);
+
+	// If rule has been changed, erase  it.
+	if (!rule_unchanged) {
+		properties.erase(rules_ptr);
+		return false;
+		} 
+	else {
+		return true;
+		}
+	}
+
+int substitute::Rules::size() {
+	return properties.size();
+}
+
+void substitute::Rules::cleanup() {
+	// Erase rules that are pointing to garbage-collected Ex expressions
+	// or rules that have possibly been changed.
+	for (auto it = properties.begin(); it != properties.end(); ) {
+		if (it->first.expired()) {
+			it = properties.erase(it);
+			} 
+		else if (it->first.lock()->state() != result_t::l_checkpointed) {
+			it = properties.erase(it);
+			}
+		else {
+			++it;
+			}
+		}
+    }

core/algorithms/substitute.hh

@@ -35,6 +35,35 @@ namespace cadabra {
 			virtual result_t apply(iterator&);
 
 			Ex_comparator comparator;
+
+			// Rules is a class for caching properties of substitution rules to avoid
+			// processing them in subsequent calls
+			class Rules {
+				public:
+					// Associate rule properties with a specific object
+					void store(Ex& rules, std::map<iterator, bool>& lhs_contains_dummies, std::map<iterator, bool>& rhs_contains_dummies);
+					// Check if rules are present
+					bool is_present(Ex& rules);
+					// Retrieve properties from rules
+					void retrieve(Ex& rules, std::map<iterator, bool>& lhs_contains_dummies, std::map<iterator, bool>& rhs_contains_dummies);
+					// Count number of rules
+					int size();
+					// Eliminate rules that are expired
+					void cleanup();
+
+				private:
+					// Map storing weak pointers to `Ex` and pairs of lhs/rhs maps as values
+					std::map<std::weak_ptr<Ex>, 
+								std::pair<std::map<iterator, bool>, std::map<iterator, bool>>, 
+								std::owner_less<std::weak_ptr<Ex>>> properties;
+
+					// Initial size threshold to trigger cleanup_rules
+					unsigned int cleanup_threshold = 100;
+					// Store max size of the rules list to avoid it getting out of hand
+					unsigned int max_size = 1000;
+				};
+
+
 		private:
 			Ex&        args;
 
@@ -48,4 +77,9 @@ namespace cadabra {
 			bool            partial;
 		};
 
+	/* Global instance of substitute::rules */
+	extern substitute::Rules replacement_rules;
+
 	}
+
+

core/pythoncdb/py_algorithms.cc

@@ -130,7 +130,7 @@ namespace cadabra {
 		def_algo<factor_in, Ex>(m, "factor_in", true, false, 0, py::arg("factors"));
 		def_algo<factor_out, Ex, bool>(m, "factor_out", true, false, 0, py::arg("factors"), py::arg("right") = false);
 		def_algo<fierz, Ex>(m, "fierz", true, false, 0, py::arg("spinors"));
-		def_algo<substitute, Ex, bool>(m, "substitute", true, false, 0, py::arg("rules"), py::arg("partial") = true);
+		def_algo<substitute, Ex&, bool>(m, "substitute", true, false, 0, py::arg("rules"), py::arg("partial") = true);
 		def_algo<take_match, Ex>(m, "take_match", true, false, 0, py::arg("rules"));
 		def_algo<replace_match>(m, "replace_match", false, false, 0);
 		def_algo<zoom, Ex, bool>(m, "zoom", true, false, 0, py::arg("rules"), py::arg("partial") = true);

web2/cadabra2/source/changelog.html

@@ -17,11 +17,11 @@ <h1>Change log</h1>
   even-numbered ones are released in packaged/installer form.
 </p>
 <a name="devel"></a>
-<h3>github devel branch (2.5.7)</h3>
+<h3>github devel branch (2.5.9)</h3>
 <ul>
 </ul>
 
-<a name="master"></a>
+<a name="2.5.8"></a>
 <h3>2.5.8 (released 25-Oct-2024)</h3>
 <ul>
   <li>Fix handling of spurious backslashes
@@ -45,13 +45,15 @@ <h3>2.5.8 (released 25-Oct-2024)</h3>
   <li>Generate source tarball asset with microtex included.</li>
 </ul>
 
+<a name="2.5.6"></a>
 <h3>2.5.6 (released 29-Sep-2024)</h3>
 <ul>
   <li>Windows build is working again, now using MSYS2 because that's the only system that provides binary packages of our dependencies and also still carries gtkmm-3.0. An installer is now auto-generated by the CI on every release.</li>
   <li>Updates to the AppImage build to fix issues with sympy and matplotlib.</li>
   <li>Many other small fixes and updates to ensure that all binaries can be built by the CI.</li>
 </ul>
 
+<a name="2.5.2"></a>
 <h3>2.5.2 (released 14-Jun-2024)</h3>
 <ul>
   <li>Fix simplification of <tt>pow</tt> nodes.</li>
@@ -74,6 +76,7 @@ <h3>2.5.2 (released 14-Jun-2024)</h3>
   <li>Added dynamical cell updates, e.g. to create animations.</li>
 </ul>
 
+<a name="2.4.4"></a>
 <h3>2.4.4 (released 20-Sep-2023)</h3>
 <ul>
   <li>Fix canonicalisation of rationals.</li>
@@ -82,6 +85,7 @@ <h3>2.4.4 (released 20-Sep-2023)</h3>
   <li>Prevent dummy indices from leaking out of <tt>\pow</tt> nodes.</li>
 </ul>
 
+<a name="2.4.2"></a>
 <h3>2.4.2 (released 22-Oct-2022)</h3>
 <ul>
   <li>Add missing trig functions for display and passthrough to sympy (Oscar).</li>
@@ -97,6 +101,7 @@ <h3>2.4.2 (released 22-Oct-2022)</h3>
   <li>Fix compatibility issues with Python-3.11.</li>
 </ul>
 
+<a name="2.4.0"></a>
 <h3>2.4.0 (released 25-Aug-2022)</h3>
 <ul>
   <li>Version to accompany the publication of

web2/cadabra2/source/index.html

@@ -26,6 +26,10 @@ <h1>News</h1>
   </div>
 </p>
 <dl>
+  <dt>25 October 2024</dt>
+  <dd>Release of 2.5.8 (<a href="changelog.html#2.5.8">changes</a>).
+    Mostly bug fixes and some improvements to the notebook interface.
+  </dd>
   <dt>24 September 2024</dt>
   <dd>Windows installer now available again, get it from the
     <a href="https://github.com/kpeeters/cadabra2/releases/latest">github release page</a>.