Merge branch 'cbehan-trace_op'

Author: kpeeters

core/Algorithm.cc

@@ -34,6 +34,7 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #include "properties/Indices.hh"
 #include "properties/Coordinate.hh"
 #include "properties/Symbol.hh"
+#include "properties/Trace.hh"
 #include "properties/DependsBase.hh"
 
 #include <sstream>
@@ -319,7 +320,8 @@ void Algorithm::propagate_zeroes(post_order_iterator& it, const iterator& topnod
 		return;
 
 	const Derivative *der=kernel.properties.get<Derivative>(walk);
-	if(*walk->name=="\\prod" || der) {
+	const Trace *trace=kernel.properties.get<Trace>(walk);
+	if(*walk->name=="\\prod" || der || trace) {
 		if(der && it->is_index()) return;
 		walk->multiplier=rat_set.insert(0).first;
 		it=walk;

core/Compare.cc

@@ -1158,6 +1158,23 @@ namespace cadabra {
 		return sign;
 		}
 
+	int Ex_comparator::can_swap_components(Ex::iterator one, Ex::iterator two, match_t subtree_comparison)
+		{
+		int tmp;
+		auto impi=properties.get_with_pattern<ImplicitIndex>(one, tmp);
+		if(impi.first) {
+			if(impi.first->explicit_form.size()>0) {
+				one=impi.first->explicit_form.begin();
+				}
+			}
+		impi=properties.get_with_pattern<ImplicitIndex>(two, tmp);
+		if(impi.first) {
+			if(impi.first->explicit_form.size()>0) {
+				two=impi.first->explicit_form.begin();
+				}
+			}
+		return can_swap(one, two, subtree_comparison, true);
+		}
 
 	// Determine the sign required to move the last factor in 'factors' to
 	// the right of the first.  If moving left, do not count signs from

core/Compare.hh

@@ -313,6 +313,14 @@ namespace cadabra {
 			int  can_swap(Ex::iterator one, Ex::iterator two, match_t subtree_comparison,
 			              bool ignore_implicit_indices=false);
 
+			/// Wrapper for can_swap which is meant for objects that have implicit
+			/// indices. This checks whether a single component of A commutes or
+			/// anticommutes with a single component of B, saying nothing about
+			/// whether A and B commute under matrix multiplication.
+
+			int  can_swap_components(Ex::iterator one, Ex::iterator two,
+						 match_t subtree_comparison);
+
 			/// Determine whether object 'one' and 'two' can be moved next
 			/// to each other by moving either one or the other: if fix_one==true
 			/// the first node is kept fixed, otherwise the second node is kept fixed.

core/Props.hh

@@ -503,7 +503,7 @@ namespace cadabra {
 				if(ret1) { // property of the right type found for object 1
 					property_map_t::const_iterator walk2=pit2.first;
 					while(walk2!=pit2.second) {
-						if((*walk2).second.first->match(*this, it2, ignore_parent_rel)) { // match for object 1 found
+						if((*walk2).second.first->match(*this, it2, ignore_parent_rel)) { // match for object 2 found
 							ret2=dynamic_cast<const T *>((*walk2).second.second);
 							if(ret2) { // property of the right type found for object 2
 								if(ret1==ret2 && walk1!=walk2) { // accept if properties are the same and patterns are not

core/algorithms/canonicalise.cc

@@ -4,6 +4,7 @@
 #include "Functional.hh"
 #include "algorithms/canonicalise.hh"
 #include "modules/xperm_new.h"
+#include "properties/Trace.hh"
 #include "properties/Traceless.hh"
 #include "properties/Diagonal.hh"
 #include "properties/Derivative.hh"
@@ -50,18 +51,53 @@ bool canonicalise::remove_traceless_traces(iterator& it)
 	while(facit!=tr.end(it)) {
 		const Traceless *trl=kernel.properties.get_composite<Traceless>(facit);
 		if(trl) {
+			unsigned int ihits=0;
 			tree_exact_less_mod_prel_obj comp(&kernel.properties);
 			std::set<Ex, tree_exact_less_mod_prel_obj> countmap(comp);
 			index_iterator indit=begin_index(facit);
 			while(indit!=end_index(facit)) {
+				bool incremented_now=false;
+				auto ind=kernel.properties.get<Indices>(indit, true);
+				if(ind) {
+					// The indexs need to be in the set for which the object is
+					// traceless (if specified, otherwise accept all).
+					if(ind->set_name==trl->index_set_name || trl->index_set_name=="") {
+						incremented_now=true;
+						++ihits;
+						}
+					}
+				else incremented_now=true;
+				// Having no name is treated as having the right name
 				if(countmap.find(Ex(indit))==countmap.end()) {
 					countmap.insert(Ex(indit));
-					++indit;
 					}
-				else {
+				else if(incremented_now) {
 					zero(it->multiplier);
 					return true;
 					}
+				++indit;
+				}
+			iterator parent=it;
+			if (tr.number_of_children(it)==1 && !tr.is_head(it)) parent=tr.parent(it);
+			const Trace *trace=kernel.properties.get<Trace>(parent);
+			if(trace) {
+				int tmp;
+				auto impi=kernel.properties.get_with_pattern<ImplicitIndex>(facit, tmp);
+				if(impi.first->explicit_form.size()>0) {
+					// Does the explicit form have two more indices of the right type?
+					Ex::iterator eform=impi.first->explicit_form.begin();
+					unsigned int ehits=0;
+					indit=begin_index(eform);
+					while(indit!=end_index(eform)) {
+						auto ind=kernel.properties.get<Indices>(indit, true);
+						if(ind->set_name==trl->index_set_name && ind->set_name==trace->index_set_name) ++ehits;
+						if(ehits - ihits > 1) {
+							zero(it->multiplier);
+							return true;
+							}
+						++indit;
+						}
+					}
 				}
 			}
 		++facit;

core/algorithms/combine.cc

@@ -5,8 +5,8 @@
 
 using namespace cadabra;
 
-combine::combine(const Kernel& k, Ex& e)
-	: Algorithm(k, e)
+combine::combine(const Kernel& k, Ex& e, Ex& t)
+	: Algorithm(k, e), trace_op(t)
 	{
 	}
 
@@ -36,7 +36,7 @@ Algorithm::result_t combine::apply(iterator& it)
 
 	while(ind_dummy.begin()!=ind_dummy.end()) {
 		bool found=false;
-		index_map_t::iterator start=ind_dummy.begin();
+		index_map_t::iterator start=ind_dummy.begin(), backup;
 		while(!found && start!=ind_dummy.end()) {
 			iterator parent=tr.parent(start->second);
 			sibling_iterator ch=tr.begin(parent), last_part;
@@ -46,11 +46,17 @@ Algorithm::result_t combine::apply(iterator& it)
 				++ch;
 				}
 			if(last_part==start->second) {
-				// We are on a rightmost contracted index
-				found=true;
+				++last_part;
+				if(last_part==tr.end(parent)) {
+					// Dummy index with nothing to the right is preferred
+					found=true;
+					}
+				else backup=start;
 				}
-			else ++start;
+			if(!found) ++start;
 			}
+		// As a backup, we use a dummy index with only non-dummies to the right
+		if(!found) start=backup;
 		bool paired=true;
 		while(paired && start!=ind_dummy.end()) {
 			iterator parent=tr.parent(start->second);
@@ -95,6 +101,7 @@ Algorithm::result_t combine::apply(iterator& it)
 			}
 		}
 
+	std::string trace_start="";
 	std::vector<Ex::iterator>::iterator dums1=dummies.begin(), dums2;
 	dums2=dums1;
 	++dums2;
@@ -133,6 +140,38 @@ Algorithm::result_t combine::apply(iterator& it)
 			iterator brackprod=tr.append_child(outerbrack, str_node("\\prod"));
 			iterator parn1=tr.parent(*dums1);
 			iterator parn2=tr.parent(*dums2);
+
+			// count how many sign changes stand between the two objects
+			int sign=1;
+			unsigned int hits=0;
+			Ex_comparator compare(kernel.properties);
+			sib=tr.begin(it);
+			while(hits<2) {
+				if(hits==1 && sib!=parn2) {
+					// pass arguments manually as can_swap() does not check them
+					bool isbrack=*(sib->name)=="\\indexbracket";
+					if(isbrack && isbrack2) {
+						auto es=compare.equal_subtree(tr.begin(parn2), tr.begin(sib));
+						sign*=compare.can_swap_components(tr.begin(parn2), tr.begin(sib), es);
+						}
+					else if(isbrack && !isbrack2) {
+						auto es=compare.equal_subtree(parn2, tr.begin(sib));
+						sign*=compare.can_swap_components(parn2, tr.begin(sib), es);
+						}
+					else if(!isbrack && isbrack2) {
+						auto es=compare.equal_subtree(tr.begin(parn2), sib);
+						sign*=compare.can_swap_components(tr.begin(parn2), sib, es);
+						}
+					else {
+						auto es=compare.equal_subtree(parn2, sib);
+						sign*=compare.can_swap_components(parn2, sib, es);
+						}
+					}
+				if(sib==parn1 || sib==parn2) ++hits;
+				++sib;
+				}
+			if(sign==-1) flip_sign(brackprod->multiplier);
+
 			// remove the dummy index from these two objects, and move
 			// other (dummy or not) indices to the outer indexbracket.
 			sibling_iterator ind1=tr.begin(tr.parent(*dums1));
@@ -197,7 +236,22 @@ Algorithm::result_t combine::apply(iterator& it)
 		if(consecutive) {
 			++dums1;
 			++dums2;
+			if(dums2!=dummies.end() && trace_op.size()>0) {
+				if(*(*dums2)->name==trace_start) {
+					iterator parn=tr.parent(*dums2);
+					iterator trace=tr.insert(parn, str_node(trace_op.begin()->name));
+					sibling_iterator nxt=tr.begin(parn);
+					++nxt;
+					++dums1;
+					++dums2;
+					tr.reparent(trace, tr.begin(parn), nxt);
+					multiply(trace->multiplier, *parn->multiplier);
+					tr.erase(parn);
+					trace_start="";
+					}
+				}
 			}
+		else trace_start=*(*dums1)->name;
 		++dums1;
 		++dums2;
 		}

core/algorithms/combine.hh

@@ -7,14 +7,15 @@ namespace cadabra {
 
 	class combine : public Algorithm {
 		public:
-			combine(const Kernel&, Ex&);
+			combine(const Kernel&, Ex&, Ex&);
 
 			virtual bool     can_apply(iterator) override;
 			virtual result_t apply(iterator&) override;
 
 		private:
 			typedef std::map<nset_t::iterator, iterator> indexlocmap_t;
 
+			Ex trace_op;
 			indexlocmap_t iloc;
 		};
 

core/algorithms/sort_product.cc

@@ -1,6 +1,7 @@
 
 #include "Cleanup.hh"
 #include "algorithms/sort_product.hh"
+#include "properties/Trace.hh"
 
 using namespace cadabra;
 
@@ -88,6 +89,61 @@ Algorithm::result_t sort_product::apply(iterator& st)
 			}
 		}
 
+	if(ret==result_t::l_no_action) {
+		bool found=false;
+		bool ascend=true;
+		iterator parn=st;
+		// It is fine to cycle when we have a sum of a sum inside a trace
+		while(ascend && !tr.is_head(parn)) {
+			parn=tr.parent(parn);
+			const Trace *trace=kernel.properties.get<Trace>(parn);
+			if(trace) {
+				ascend=false;
+				found=true;
+				}
+			else if(*parn->name=="\\indexbracket") {
+				ascend=false;
+				if(number_of_indices(parn)==2) {
+					index_iterator first=begin_index(parn);
+					index_iterator last=first;
+					++last;
+					if(*first->name==*last->name) found=true;
+					}
+				}
+			else if(*parn->name!="\\sum") {
+				ascend=false;
+				}
+			}
+		if(found) {
+			one=tr.begin(st);
+			two=one;
+			++two;
+			while(two!=tr.end(st)) {
+				compare.clear();
+				auto es=compare.equal_subtree(one, two);
+				if(compare.should_swap(one, es)) one=two;
+				++two;
+				}
+			// We have found the element that should go at the front of the trace
+			Ex::sibling_iterator front=one;
+			while(tr.begin(st)!=front) {
+				one=tr.begin(st);
+				two=one;
+				++two;
+				while(two!=tr.end(st)) {
+					compare.clear();
+					auto es=compare.equal_subtree(one, two);
+					int sign=compare.can_swap_components(one, two, es);
+					if(sign==-1) flip_sign(st->multiplier);
+					tr.swap(one);
+					++two;
+					++two;
+					}
+				}
+			ret=result_t::l_applied;
+			}
+		}
+
 	if(cleanup)
 		cleanup_dispatch(kernel, tr, st);
 

core/algorithms/untrace.cc

@@ -75,8 +75,16 @@ Algorithm::result_t untrace::apply(iterator& trloc)
 			}
 
 		if(move_out) {
+			int sign=1;
+			sibling_iterator st2=tr.begin(prodloc);
+			Ex_comparator compare(kernel.properties);
+			while(st2!=st) {
+				auto es=compare.equal_subtree(st, st2);
+				sign*=compare.can_swap_components(st, st2, es);
+				++st2;
+				}
 			tr.move_before(trloc, st);
-			//multiply(it->multiplier, sign);
+			multiply(trloc->multiplier, sign);
 			}
 
 		st=nxt;

core/properties/Traceless.cc

@@ -7,3 +7,11 @@ std::string Traceless::name() const
 	{
 	return "Traceless";
 	}
+
+bool Traceless::parse(Kernel&, keyval_t& keyvals)
+	{
+	keyval_t::const_iterator kv=keyvals.find("indices");
+	if(kv!=keyvals.end())
+		index_set_name=*(kv->second->name);
+	return true;
+	}