Merge branch 'fix/evaluate_with_partial'

Author: kpeeters

core/DisplayTeX.cc

@@ -430,8 +430,8 @@ void DisplayTeX::print_fraclike(std::ostream& str, Ex::iterator it)
 	if(mult * (*it->multiplier)!=1) {
 		print_multiplier(str, it, mult);
 		}
-
-	dispatch(str, num);
+	if(num->is_rational()==false || (mult * (*it->multiplier))==1)
+		dispatch(str, num);
 	str << "}{";
 	dispatch(str, den);
 	str << "}";

core/DisplayTerminal.cc

@@ -352,7 +352,8 @@ void DisplayTerminal::print_fraclike(std::ostream& str, Ex::iterator it)
 //	if(needs_brackets(num))
 //		str << "(";
 
-	dispatch(str, num);
+	if(num->is_rational()==false || (*it->multiplier)==1)
+		dispatch(str, num);
 
 //	if(needs_brackets(num))
 //		str << ")";

core/Permutations.hh

@@ -24,10 +24,11 @@ class Perm {
 		std::vector<int> perm;
 
 		/// Find the permutation that takes [start1, end1> to [start2, end2>.
+		/// This will be available in 'perm' afterwards.
 		template<class iterator>
 		void find(iterator start1, iterator end1, iterator start2, iterator end2);
 
-		/// Apply the permutation on the range. 
+		/// Apply the permutation 'perm' on the given range.
 		template<class iterator>
 		void apply(iterator start1, iterator end1);
 
@@ -81,7 +82,11 @@ void Perm::apply(iterator start, iterator end)
 		}
 
 	// std::cerr << orig.size() << ", " << perm.size() << std::endl;
-	assert(orig.size()==perm.size());
+	if(orig.size()!=perm.size()) {
+		std::cerr << "Perm::apply: orig.size()=" << orig.size() << ", "
+					 << "perm.size()=" << perm.size() << std::endl;
+		assert(orig.size()==perm.size());
+		}
 
 	for(unsigned int i=0; i<orig.size(); ++i) {
 		*start=orig[perm[i]];

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

@@ -364,18 +364,24 @@ void evaluate::merge_components(iterator it1, iterator it2)
 	// we can be assured that the free indices match; they just may not be
 	// in the same order).
 
+	#ifdef DEBUG
+	std::cerr << "merge_components on " << Ex(it1) << " and " << Ex(it2) << std::endl;
+	#endif
+
 	assert(*it1->name=="\\components");
 	assert(*it2->name=="\\components");
 	sibling_iterator sib1=tr.end(it1);
 	--sib1;
 	sibling_iterator sib2=tr.end(it2);
 	--sib2;
+	assert(*sib1->name=="\\comma");
+	assert(*sib2->name=="\\comma");	
 
-	// We cannot directly compare the lhs of this equals node with the lhs
-	// of the equals node of the other components node, because the index
-	// order on the two components nodes may be different. We first
-	// have to ensure that the orders are the same (but only, of course)
-	// if we have anything to permutate in the first place.
+	// We cannot directly compare the lhs of the equals nodes of it1
+	// with the lhs of the equals node of it2, because the index order
+	// on the two components nodes may be different. We first have to
+	// ensure that the orders are the same (but only, of course) if we
+	// have anything to permutate in the first place.
 
 	if(*tr.begin(it1)->name!="\\comma") {
 		// Look at all indices on the two components nodes. Find
@@ -384,14 +390,20 @@ void evaluate::merge_components(iterator it1, iterator it2)
 		Perm perm;
 		perm.find(tr.begin(it2), sib2, tr.begin(it1), sib1);
 
-		//perm.apply(tr.begin(it2), sib2);
-		//std::cerr << "after permutation " << Ex(tr) << std::endl;
-		
+		// For each \equals node in the it2 comma node, permute
+		// the values so they agree with the index order on it1.
 		cadabra::do_list(tr, sib2, [&](Ex::iterator nd) {
-				auto lhs2 = tr.begin(nd);
-				perm.apply(tr.begin(lhs2), tr.end(lhs2));
+				// nd is an \equals node.
+				assert(*nd->name=="\\equals");
+				auto comma = tr.begin(nd);
+				assert(*comma->name=="\\comma");
+				perm.apply(tr.begin(comma), tr.end(comma));
 				return true;
 				});
+
+      #ifdef DEBUG
+		std::cerr << "permutations done" << std::endl;
+		#endif
 		}
 
 	// Now all index orders match and we can simply compare index value sets.
@@ -527,8 +539,10 @@ Ex::iterator evaluate::handle_derivative(iterator it)
 	size_t ni=number_of_direct_indices(it);
 
 	cadabra::do_list(tr, ivalues, [&](Ex::iterator iv) {
-			// std::cerr << "====" << std::endl;
-			// std::cerr << Ex(iv) << std::endl;
+			#ifdef DEBUG
+			std::cerr << "====" << std::endl;
+			std::cerr << Ex(iv) << std::endl;
+			#endif
 			// For each internal dummy set, keep track of the
 			// position in the permutation array where we generate
 			// its value.
@@ -539,9 +553,11 @@ Ex::iterator evaluate::handle_derivative(iterator it)
 			auto deps=dependencies(rhs);
 
 			// If the argument does not depend on anything, all derivatives
-			// would produce zero.
-			if(deps.size()==0) 
+			// would produce zero. Remove this \equals node from the tree.
+			if(deps.size()==0) {
+				tr.erase(iv);
 				return true;
+				}
 
 			// All indices on \partial can take any of the values of the
 			// dependencies, EXCEPT when the index is a dummy index. In
@@ -556,7 +572,9 @@ Ex::iterator evaluate::handle_derivative(iterator it)
 
 			combin::combinations<Ex> cb;
 			for(auto& obj: deps) {
-				// std::cerr << "dep " << obj << std::endl;
+				#ifdef DEBUG
+				std::cerr << "dep " << obj << std::endl;
+				#endif
 				cb.original.push_back(obj);
 				}
 			cb.multiple_pick=true;
@@ -596,18 +614,24 @@ Ex::iterator evaluate::handle_derivative(iterator it)
 			// derivative, create an entry in the \components node.
 
 			for(unsigned int i=0; i<cb.size() || cb.size()==0; ++i) {
-				// std::cerr << "Index combination " << i << std::endl;
+				#ifdef DEBUG
+				std::cerr << "Index combination " << i << std::endl;
+				#endif
 				Ex eqcopy(iv); 
 				auto lhs=eqcopy.begin(eqcopy.begin());
 				assert(*lhs->name=="\\comma");
 
 				if(cb.size()>0) {
+					#ifdef DEBUG
+					std::cerr << "Copying values of derivative indices" << std::endl;
+					#endif
 					// Setup the index values; simply copy from the cb array, but only
 					// if the indices are not internal dummy.
 					for(size_t j=0; j<cb[i].size(); ++j) {
 						auto fd = ind_dummy.find(Ex(tr.child(it, j)));
-						if(fd==ind_dummy.end()) 
+						if(fd==ind_dummy.end()) {
 							eqcopy.append_child(iterator(lhs), cb[i][j].begin() );
+							}
 						}
 					}
 				auto rhs=lhs;
@@ -672,14 +696,16 @@ Ex::iterator evaluate::handle_derivative(iterator it)
 
 				if(cb.size()==0) break;
 				}
+
+			// Erase the original \equals entry (we generated a full replacement above).
 			tr.erase(iv);
 			return true;
 			});
 
    one(it->multiplier);
 	// std::cerr << "now " << Ex(it) << std::endl;
 
-	
+
 	// Now move the free (but not the internal dummy!) partial indices
    //	to the components node, and then unwrap the partial node.
 
@@ -707,11 +733,15 @@ Ex::iterator evaluate::handle_derivative(iterator it)
 			++se;
 		}
 
-		// std::cerr << "after index move " << Ex(it) << std::endl;
+	#ifdef DEBUG
+	std::cerr << "after index move " << Ex(it) << std::endl;
+	#endif
 
 	merge_component_children(it);
 
-	// std::cerr << "after merge " << Ex(it) << std::endl;
+	#ifdef DEBUG
+	std::cerr << "after merge " << Ex(it) << std::endl;
+	#endif
 
 	simplify_components(it);
 	// std::cerr << "then " << Ex(it) << std::endl;

doc/cadabra2_hep.tex

@@ -493,9 +493,9 @@ \subsection{Bianchi, Ricci and Schouten identities}
 
 \toprule
 \begin{cdbin}
-\tableau{#}::FilledTableau(dimension=10).
-\tableau{0,0}{1,1} \tableau{a,a}{b,b}:
-@lr_tensor!(%);
+\ftableau{#}::FilledTableau(dimension=10).
+ex:=\ftableau{0,0}{1,1} \ftableau{a,a}{b,b}:
+lr_tensor(_);
 \end{cdbin}
 \begin{cdbout}
 \ftableau{{0}{0}{a}{a},{1}{1}{b}{b}} + \ftableau{{0}{0}{a}{a},{1}{1}{b},{b}} + \ftableau{{0}{0}{a}{a},{1}{1},{b}{b}} + \ftableau{{0}{0}{a},{1}{1}{b},{a},{b}} + \ftableau{{0}{0}{a},{1}{1},{a}{b},{b}} + \ftableau{{0}{0},{1}{1},{a}{a},{b}{b}};

tests/components.cdb

@@ -418,6 +418,7 @@ def test22():
 test22()
 
 def test23():
+    __cdbkernel__=create_scope()
     {r,t}::Coordinate;        
     {a,b,c}::Indices(values={r,t}, position=fixed);
     \partial{#}::PartialDerivative;
@@ -453,18 +454,42 @@ test23()
 # rl:= \partial_{r}{ A_{t} } = 1;
 # evaluate(ex, rl);
 
-{\theta, \varphi}::Coordinate;
-{\alpha, \beta, \gamma, \delta, \rho, \sigma, \mu, \nu, \lambda}::Indices(values={\varphi, \theta}, position=fixed);
-\partial{#}::PartialDerivative;
-g_{\alpha\beta}::Metric.
-g^{\alpha\beta}::InverseMetric.
+def test24():
+    __cdbkernel__=create_scope()
+    {\theta, \varphi}::Coordinate;
+    {\alpha, \beta, \gamma, \delta, \rho, \sigma, \mu, \nu, \lambda}::Indices(values={\varphi, \theta}, position=fixed);
+    \partial{#}::PartialDerivative;
+    g_{\alpha\beta}::Metric.
+    g^{\alpha\beta}::InverseMetric.
+
+    sphe:={ g_{\theta\theta} = r**2,
+        g_{\varphi\varphi} = r**2 \sin(\theta)**2 }.
+    complete(sphe, $g^{\alpha\beta}$)
+
+    ch:= \Gamma^{\alpha}_{\mu\nu} = 1/2 g^{\alpha\beta} ( \partial_{\nu}{g_{\beta\mu}}
+                                                        +\partial_{\mu}{g_{\beta\nu}}
+                                                        -\partial_{\beta}{g_{\mu\nu}} );
 
-sphe:={ g_{\theta\theta} = r**2,
-         g_{\varphi\varphi} = r**2 \sin(\theta)**2 }.
-# complete(sphe, $g^{\alpha\beta}$);
+    evaluate(ch, sphe, rhsonly=True)
 
-#ch:= \Gamma^{\alpha}_{\mu\nu} = 1/2 g^{\alpha\beta} ( \partial_{\nu}{g_{\beta\mu}}
-#                                                      +\partial_{\mu}{g_{\beta\nu}}
-#                                                     -\partial_{\beta}{g_{\mu\nu}} );
+    rm:= R^{\rho}_{\sigma\mu\nu} = +\partial_{\mu}{\Gamma^{\rho}_{\sigma\nu}}
+           -\partial_{\nu}{\Gamma^{\rho}_{\sigma\mu}}
+           +\Gamma^{\rho}_{\beta\mu} \Gamma^{\beta}_{\sigma\nu}
+           -\Gamma^{\rho}_{\beta\nu} \Gamma^{\beta}_{\sigma\mu};
 
-# evaluate(ch, sphe, rhsonly=True);
+    substitute(rm, ch)
+    evaluate(rm, sphe, rhsonly=True)
+
+    ricci:= R_{\sigma\nu} = R^{\rho}_{\sigma\rho\nu};
+    substitute(ricci, rm)
+    evaluate(ricci, sphe, rhsonly=True)
+
+    R:= R = R_{\sigma\nu} g^{\sigma\nu};
+    substitute(R, ricci)
+    evaluate(R, sphe, rhsonly=True)
+    val = rhs(R)
+    tst:= @(val) - 2/r**2;
+    assert(tst==0)
+    print("Test 24 passed")
+        
+test24()