Finding bug in map_sympy with tensor prefactors

Author: kpeeters

core/SympyCdb.cc

@@ -11,7 +11,7 @@
 
 using namespace cadabra;
 
-//#define DEBUG
+#define DEBUG
 
 Ex::iterator sympy::apply(const Kernel& kernel, Ex& ex, Ex::iterator& it, const std::vector<std::string>& wrap, const std::string& args, 
 								  const std::string& method)

core/algorithms/evaluate.cc

@@ -159,7 +159,7 @@ Ex::iterator evaluate::handle_sum(iterator it)
 	auto sib2=sib1;
 	++sib2;
 	while(sib2!=tr.end(it)) {
-		std::cerr << "merging components " << Ex(sib1) << " and " << Ex(sib2) << std::endl;
+//		std::cerr << "merging components " << Ex(sib1) << " and " << Ex(sib2) << std::endl;
 		merge_components(sib1, sib2);
 		sib2=tr.erase(sib2);
 		}

core/algorithms/map_sympy.cc

@@ -35,6 +35,9 @@ bool map_sympy::can_apply(iterator st)
    // In a product, it is still possible that there is a sub-product which
 	// contains no indices.
 	if(*st->name=="\\prod") {
+		// Find the factors in the product which have a proper index on them. Do this by
+		// starting at the index, and if it is not coordinate or symbol, then go up until we
+		// reach the first child level of the product.
 		for(auto& ind: ind_free) {
 			const Coordinate *cdn=kernel.properties.get_composite<Coordinate>(ind.second, true);
 			const Symbol     *smb=kernel.properties.get_composite<Symbol>(ind.second, true);
@@ -73,13 +76,20 @@ Algorithm::result_t map_sympy::apply(iterator& it)
 	wrap.push_back(head_);
 
 	if(left.size()>0) {
-		std::cerr << "Sub-product with " << left.size() << " non-index carrying factors" << std::endl;
 		Ex prod("\\prod");
 		for(auto& fac: left)
-			prod.append_child(fac);
+			prod.append_child(prod.begin(), fac);
 		auto top=prod.begin();
 		sympy::apply(kernel, prod, top, wrap, "", "");
-		for(auto& kl: index_factors)
+		// Now remove the non-index carrying factors and replace with
+		// the factors of 'prod' just simplified.
+		sibling_iterator ps=prod.begin(top);
+		while(ps!=prod.end(top)) {
+			tr.insert_subtree(*left.begin(), ps);
+			++ps;
+			}
+		std::cerr << "Before erasing " << Ex(it) << std::endl;
+		for(auto& kl: left)
 			tr.erase(kl);
 
 		return result_t::l_no_action;