Add Determinant and Trace properties and related logic in complete.

Author: kpeeters

core/CMakeLists.txt

@@ -103,6 +103,7 @@ set(PROPERTY_SRC_FILES
 	properties/Depends.cc
 	properties/DependsInherit.cc
 	properties/Derivative.cc
+	properties/Determinant.cc
 	properties/Diagonal.cc
 	properties/DifferentialForm.cc
 	properties/DiracBar.cc
@@ -136,6 +137,7 @@ set(PROPERTY_SRC_FILES
 	properties/Tableau.cc
 	properties/TableauBase.cc
 	properties/TableauSymmetry.cc
+	properties/Trace.cc	
 	properties/Traceless.cc
 	properties/Vielbein.cc
 	properties/Weight.cc

core/PythonCdb.cc

@@ -38,6 +38,7 @@
 #include "properties/Depends.hh"
 #include "properties/DependsInherit.hh"
 #include "properties/Derivative.hh"
+#include "properties/Determinant.hh"
 #include "properties/DifferentialForm.hh"
 #include "properties/DiracBar.hh"
 #include "properties/GammaMatrix.hh"
@@ -73,6 +74,7 @@
 #include "properties/Symmetric.hh"
 #include "properties/Tableau.hh"
 #include "properties/TableauSymmetry.hh"
+#include "properties/Trace.hh"
 #include "properties/Traceless.hh"
 #include "properties/Vielbein.hh"
 #include "properties/Weight.hh"
@@ -1569,6 +1571,7 @@ PYBIND11_MODULE(cadabra2, m)
 	def_prop<DAntiSymmetric>(m);
 	def_prop<Depends>(m);
 	def_prop<Derivative>(m);
+	def_prop<Determinant>(m);	
 	def_prop<Diagonal>(m);
 	def_prop<DifferentialForm>(m);
 	def_prop<Distributable>(m);
@@ -1601,6 +1604,7 @@ PYBIND11_MODULE(cadabra2, m)
 	def_prop<Symmetric>(m);
 	def_prop<Tableau>(m);
 	def_prop<TableauSymmetry>(m);
+	def_prop<Trace>(m);	
 	def_prop<Traceless>(m);
 	def_prop<Vielbein>(m);
 	def_prop<InverseVielbein>(m);

core/SympyCdb.cc

@@ -12,7 +12,7 @@
 
 using namespace cadabra;
 
-//#define DEBUG 1
+// #define DEBUG 1
 
 Ex::iterator sympy::apply(const Kernel& kernel, Ex& ex, Ex::iterator& it, const std::vector<std::string>& wrap, std::vector<std::string> args, 
 								  const std::string& method)
@@ -93,7 +93,7 @@ Ex::iterator sympy::apply(const Kernel& kernel, Ex& ex, Ex::iterator& it, const
 	return it;
 	}
 
-Ex sympy::invert_matrix(const Kernel& kernel, Ex& ex, Ex& rules)
+Ex sympy::fill_matrix(const Kernel& kernel, Ex& ex, Ex& rules)
 	{
 	// check that object has two children only.
 	if(ex.number_of_children(ex.begin())!=2) {
@@ -103,12 +103,6 @@ Ex sympy::invert_matrix(const Kernel& kernel, Ex& ex, Ex& rules)
 	Ex::iterator ind1=ex.child(ex.begin(), 0);
 	Ex::iterator ind2=ex.child(ex.begin(), 1);
 
-	str_node::parent_rel_t prel;
-	if(ind1->fl.parent_rel==str_node::p_sub)   prel=str_node::p_super;
-	if(ind1->fl.parent_rel==str_node::p_super) prel=str_node::p_sub;
-
-	Ex ret;
-
 	// Get Indices property and from there Coordinates.
 
 	const Indices *prop1 = kernel.properties.get<Indices>(ind1);
@@ -130,10 +124,8 @@ Ex sympy::invert_matrix(const Kernel& kernel, Ex& ex, Ex& rules)
 			Ex c(ex.begin());
 			Ex::iterator cit1=c.child(c.begin(), 0);
 			Ex::iterator cit2=c.child(c.begin(), 1);
-			cit1=c.replace(cit1, prop1->values[c1].begin());
-			cit2=c.replace(cit2, prop1->values[c2].begin());
-			cit1->fl.parent_rel=prel;
-			cit2->fl.parent_rel=prel;
+			cit1=c.replace_index(cit1, prop1->values[c1].begin(), true);
+			cit2=c.replace_index(cit2, prop1->values[c2].begin(), true);
 
 			Ex::iterator cit=c.begin();
 			substitute subs(kernel, c, rules);
@@ -146,38 +138,80 @@ Ex sympy::invert_matrix(const Kernel& kernel, Ex& ex, Ex& rules)
 				}
 			}
 		}
+
+	return matrix;
+	}
+
+void sympy::invert_matrix(const Kernel& kernel, Ex& ex, Ex& rules, const Ex& tocompute)
+	{
+	if(ex.number_of_children(ex.begin())!=2) {
+		throw ConsistencyException("Object should have exactly two indices.");
+		}
 
+	auto matrix = fill_matrix(kernel, ex, rules);
+
 	auto top=matrix.begin();
 	std::vector<std::string> wrap;
 	sympy::apply(kernel, matrix, top, wrap, std::vector<std::string>(), ".inv()");
 	//matrix.print_recursive_treeform(std::cerr, top);
 
+	Ex::iterator ind1=ex.child(ex.begin(), 0);
+	Ex::iterator ind2=ex.child(ex.begin(), 1);
+	const Indices *prop1 = kernel.properties.get<Indices>(ind1);
+	const Indices *prop2 = kernel.properties.get<Indices>(ind2);
+
 	Ex::iterator ruleslist=rules.begin();
 
 	// Now we need to iterate over the components again and construct sparse rules.
-	cols=matrix.begin(matrix.begin()); // outer comma
+	auto cols=matrix.begin(matrix.begin()); // outer comma
 	auto row=matrix.begin(cols); // first inner comma
 	for(unsigned c1=0; c1<prop1->values.size(); ++c1) {
 		auto el =matrix.begin(row);  // first element of first inner comma
 		for(unsigned c2=0; c2<prop2->values.size(); ++c2) {
 			if(el->is_zero()==false) {
 				Ex rule("\\equals");
-				auto rit  = rule.append_child(rule.begin(), ex.begin());
+				auto rit  = rule.append_child(rule.begin(), tocompute.begin());
 				auto cvit = rule.append_child(rule.begin(), Ex::iterator(el));
 				auto i = rule.begin(rit);
 				//std::cerr << c1 << ", " << c2 << std::endl;
-				i = rule.replace(i, prop1->values[c1].begin());
-				i->fl.parent_rel=ind1->fl.parent_rel;
+				i = rule.replace_index(i, prop1->values[c1].begin(), true);
+//				i->fl.parent_rel=ind1->fl.parent_rel;
 				++i;
-				i = rule.replace(i, prop1->values[c2].begin());
-				i->fl.parent_rel=ind1->fl.parent_rel;
+				i = rule.replace_index(i, prop1->values[c2].begin(), true);
+//				i->fl.parent_rel=ind1->fl.parent_rel;
 				rules.append_child(ruleslist, rule.begin());
 				//rule.print_recursive_treeform(std::cerr, rule.begin());
 				}
 			++el;
 			}
 		++row;
 		}
+	}
+
+void sympy::determinant(const Kernel& kernel, Ex& ex, Ex& rules, const Ex& tocompute)
+	{
+	auto matrix = fill_matrix(kernel, ex, rules);
+
+	auto top=matrix.begin();
+	std::vector<std::string> wrap;
+	sympy::apply(kernel, matrix, top, wrap, std::vector<std::string>(), ".det()");
+
+	Ex rule("\\equals");
+	rule.append_child(rule.begin(), tocompute.begin());
+	rule.append_child(rule.begin(), matrix.begin());
+	rules.append_child(rules.begin(), rule.begin());
+	}
+
+void sympy::trace(const Kernel& kernel, Ex& ex, Ex& rules, const Ex& tocompute)
+	{
+	auto matrix = fill_matrix(kernel, ex, rules);
+
+	auto top=matrix.begin();
+	std::vector<std::string> wrap;
+	sympy::apply(kernel, matrix, top, wrap, std::vector<std::string>(), ".tr()");
 
-	return ret;
+	Ex rule("\\equals");
+	rule.append_child(rule.begin(), tocompute.begin());
+	rule.append_child(rule.begin(), matrix.begin());
+	rules.append_child(rules.begin(), rule.begin());
 	}

core/SympyCdb.hh

@@ -42,8 +42,17 @@ namespace sympy {
    /// sparse components. Will return a set of Cadabra rules for the
    /// inverse matrix.
 
-	cadabra::Ex invert_matrix(const cadabra::Kernel&, cadabra::Ex& ex, cadabra::Ex& rules);
+	void invert_matrix(const cadabra::Kernel&, cadabra::Ex& ex, cadabra::Ex& rules, const cadabra::Ex& tocompute);
 
+	/// \ingroup scalar
+	///
+	/// Use Sympy to compute the determinant of a matrix, given a set of rules determining
+	/// its sparse components. Will add the rules to the list.
+
+	void determinant(const cadabra::Kernel&, cadabra::Ex& ex, cadabra::Ex& rules, const cadabra::Ex& tocompute);
+	void trace(const cadabra::Kernel&, cadabra::Ex& ex, cadabra::Ex& rules, const cadabra::Ex& tocompute);
+
+	cadabra::Ex fill_matrix(const cadabra::Kernel&, cadabra::Ex& ex, cadabra::Ex& rules);
 //	extern Stopwatch sympy_stopwatch;
 
 };

core/algorithms/complete.cc

@@ -3,6 +3,8 @@
 #include "Functional.hh"
 #include "properties/Metric.hh"
 #include "properties/InverseMetric.hh"
+#include "properties/Determinant.hh"
+#include "properties/Trace.hh"
 #include "SympyCdb.hh"
 
 using namespace cadabra;
@@ -26,7 +28,24 @@ Algorithm::result_t complete::apply(iterator& it)
 	const InverseMetric *invmetric = kernel.properties.get<InverseMetric>(bg);
 	if(invmetric) {
 		Ex metric(bg);
-      sympy::invert_matrix(kernel, metric, tr);
+		Ex::iterator ind1=metric.child(metric.begin(), 0);
+		Ex::iterator ind2=metric.child(metric.begin(), 1);
+		ind1->flip_parent_rel();
+		ind2->flip_parent_rel();		
+		
+      sympy::invert_matrix(kernel, metric, tr, bg);
+		res = result_t::l_applied;
+		}
+	const Determinant *det = kernel.properties.get<Determinant>(bg);
+	if(det) {
+		Ex metric(det->obj);
+		sympy::determinant(kernel, metric, tr, bg);
+		res = result_t::l_applied;
+		}
+	const Trace *trace = kernel.properties.get<Trace>(bg);
+	if(trace) {
+		Ex metric(trace->obj);
+		sympy::trace(kernel, metric, tr, bg);
 		res = result_t::l_applied;
 		}
 

core/properties/Determinant.cc

@@ -0,0 +1,36 @@
+
+#include "Exceptions.hh"
+#include "properties/Determinant.hh"
+
+using namespace cadabra;
+
+Determinant::Determinant()
+	{
+	}
+
+std::string Determinant::name() const
+	{
+	return "Determinant";
+	}
+
+std::string Determinant::unnamed_argument() const
+	{
+	return "object";
+	}
+
+bool Determinant::parse(Kernel&, keyval_t& keyvals) 
+	{
+	keyval_t::const_iterator kv=keyvals.find("object");
+	if(kv!=keyvals.end())
+		obj = kv->second;
+	return true;
+	}
+
+void Determinant::validate(const Kernel&, const Ex& tr) const
+	{
+	}
+
+void Determinant::latex(std::ostream& str) const
+	{
+	str << name();
+	}

core/properties/Determinant.cnb

@@ -0,0 +1,110 @@
+{
+	"cells" : 
+	[
+		{
+			"cell_origin" : "client",
+			"cell_type" : "latex",
+			"cells" : 
+			[
+				{
+					"cell_origin" : "client",
+					"cell_type" : "latex_view",
+					"source" : "\\property{Determinant}{Declares an object to be the determinant of a matrix or metric.}\n\nIn order to declare the name associated to the determinant of a matrix or metric, give it the\n\\prop{Determinant} property, indicating the matrix or metric as an argument. This is mainly used\nto compute determinants in components, but the property can also be used by algorithms which\nonly deal with purely abstract properties, e.g.~covariant derivatives.\n\nThe following example shows how to declare a metric and its components, and then compute\nthe determinant from there."
+				}
+			],
+			"hidden" : true,
+			"source" : "\\property{Determinant}{Declares an object to be the determinant of a matrix or metric.}\n\nIn order to declare the name associated to the determinant of a matrix or metric, give it the\n\\prop{Determinant} property, indicating the matrix or metric as an argument. This is mainly used\nto compute determinants in components, but the property can also be used by algorithms which\nonly deal with purely abstract properties, e.g.~covariant derivatives.\n\nThe following example shows how to declare a metric and its components, and then compute\nthe determinant from there."
+		},
+		{
+			"cell_origin" : "client",
+			"cell_type" : "input",
+			"cells" : 
+			[
+				{
+					"cell_origin" : "server",
+					"cell_type" : "latex_view",
+					"source" : "\\begin{dmath*}{}\\text{Attached property Coordinate to~}\\left[t,~\\discretionary{}{}{} r,~\\discretionary{}{}{} \\phi\\right].\\end{dmath*}"
+				},
+				{
+					"cell_origin" : "server",
+					"cell_type" : "latex_view",
+					"source" : "\\begin{dmath*}{}\\text{Attached property Indices(position=fixed) to~}\\left[m,~\\discretionary{}{}{} n,~\\discretionary{}{}{} p,~\\discretionary{}{}{} q\\right].\\end{dmath*}"
+				},
+				{
+					"cell_origin" : "server",
+					"cell_type" : "latex_view",
+					"source" : "\\begin{dmath*}{}\\text{Attached property Metric to~}g_{m n}.\\end{dmath*}"
+				},
+				{
+					"cell_origin" : "server",
+					"cell_type" : "latex_view",
+					"cells" : 
+					[
+						{
+							"cell_origin" : "server",
+							"cell_type" : "input_form",
+							"source" : "{g_{t t} = -1, g_{r r} = 1, g_{\\phi \\phi} = (r)**2}"
+						}
+					],
+					"source" : "\\begin{dmath*}{}\\left[g_{t t} = -1,~\\discretionary{}{}{} g_{r r} = 1,~\\discretionary{}{}{} g_{\\phi \\phi} = {r}^{2}\\right]\\end{dmath*}"
+				}
+			],
+			"source" : "{t,r,\\phi}::Coordinate;\n{m,n,p,q}::Indices(position=fixed, values={t,r,\\phi});\ng_{m n}::Metric;\nrl:= g_{t t} = -1, g_{r r} = 1, g_{\\phi \\phi} = r**2;"
+		},
+		{
+			"cell_origin" : "client",
+			"cell_type" : "latex",
+			"cells" : 
+			[
+				{
+					"cell_origin" : "client",
+					"cell_type" : "latex_view",
+					"source" : "Now we declare that $g$ denotes the determinant of the metric $g_{m n}$, and use \\algo{complete} to compute it\ngiven the metric components."
+				}
+			],
+			"hidden" : true,
+			"source" : "Now we declare that $g$ denotes the determinant of the metric $g_{m n}$, and use \\algo{complete} to compute it\ngiven the metric components."
+		},
+		{
+			"cell_origin" : "client",
+			"cell_type" : "input",
+			"cells" : 
+			[
+				{
+					"cell_origin" : "server",
+					"cell_type" : "latex_view",
+					"source" : "\\begin{dmath*}{}\\text{Attached property Determinant to~}g.\\end{dmath*}"
+				}
+			],
+			"source" : "g::Determinant(g_{m n});"
+		},
+		{
+			"cell_origin" : "client",
+			"cell_type" : "input",
+			"cells" : 
+			[
+				{
+					"cell_origin" : "server",
+					"cell_type" : "latex_view",
+					"cells" : 
+					[
+						{
+							"cell_origin" : "server",
+							"cell_type" : "input_form",
+							"source" : "{g_{t t} = -1, g_{r r} = 1, g_{\\phi \\phi} = (r)**2, g = -(r)**2}"
+						}
+					],
+					"source" : "\\begin{dmath*}{}\\left[g_{t t} = -1,~\\discretionary{}{}{} g_{r r} = 1,~\\discretionary{}{}{} g_{\\phi \\phi} = {r}^{2},~\\discretionary{}{}{} g = -{r}^{2}\\right]\\end{dmath*}"
+				}
+			],
+			"source" : "complete(rl, $g$);"
+		},
+		{
+			"cell_origin" : "client",
+			"cell_type" : "input",
+			"source" : ""
+		}
+	],
+	"description" : "Cadabra JSON notebook format",
+	"version" : 1
+}