Inverse matrix (<4x4)

Author: lukas-valenta-tul

include/array.hh

@@ -1443,6 +1443,43 @@ public:
 		}
 	}
 
+	/*static Array det(const Array& a) {
+
+	}*/
+
+	//Square matrix inverse
+	static Array inv(const Array& a) {
+		switch (a.shape().size()) {
+		case 0: //scalar
+			Throw() << "Cannot inverse a scalar" << "\n";
+			break;
+		case 1: //vector
+		{
+			Throw() << "Cannot inverse a vector" << "\n";
+		}
+		case 2: //matrix
+		{
+			if (a.shape()[0] != a.shape()[1]) {
+				Throw() << "Cannot inverse non-square matrix" << "\n";
+			}
+			WrappedArray m_a(wrap_array(a));
+			switch (a.shape()[0]) {
+			case 1:
+				return unwrap_array(Eigen::Ref<Eigen::Matrix<details::ScalarWrapper,1,1>>(m_a).inverse());
+			case 2:
+				return unwrap_array(Eigen::Ref<Eigen::Matrix2<details::ScalarWrapper>>(m_a).inverse());
+			case 3:
+				return unwrap_array(Eigen::Ref<Eigen::Matrix3<details::ScalarWrapper>>(m_a).inverse());
+			case 4:
+				return unwrap_array(Eigen::Ref<Eigen::Matrix4<details::ScalarWrapper>>(m_a).inverse());
+			}
+			//return unwrap_array(m_a.inverse());
+		}
+		default:
+			Throw() << "Cannot inverse ND tensors" << "\n";
+		}
+	}
+
 
 	static Array flatten(const Array &tensor) {
 		uint n_elements = shape_size(tensor.shape());

include/grammar.impl.hh

@@ -193,6 +193,8 @@ struct grammar : qi::grammar<Iterator, ast::operand(), ascii::space_type> {
             FN("cross"  , &Array::cross)
             FN("sym"    , &Array::sym)
             FN("dev"    , &Array::dev)
+            //FN("det"    , &Array::det)
+            FN("inv"    , &Array::inv)
             ;
 
         unary_op.add

include/scalar_wrapper.hh

@@ -13,6 +13,7 @@
 #include "scalar_node.hh"
 #include <Eigen/Core>
 #include <Eigen/Geometry>
+#include <Eigen/LU>
 //#include <Eigen/Eigenvalues> //impossible
 
 namespace bparser {

test/test_parser.cc

@@ -275,6 +275,15 @@ void test_expression() {
 	BP_ASSERT(test_expr("dev([[1,2],[3,4]])", { 1-2.5, 2, 3, 4-2.5 }, {2,2}));
 	BP_ASSERT(test_expr("tr(dev([[1,2],[3,4]]))", { 0 }, {}));
 
+	//BP_ASSERT(test_expr("det()", {}, {}));
+
+	BP_ASSERT(test_expr("inv([[1,2],[3,4]])", { -2., 1., 1.5, -0.5 }, { 2,2 }));
+	BP_ASSERT(test_expr("a=[[1]];											 inv(a)", { 1 }, { 1,1 }));
+	BP_ASSERT(test_expr("a=[[1,2],[3,4]];									 a @ inv(a)", { 1,0,0,1 }, { 2,2 }));
+	BP_ASSERT(test_expr("a=[[1,2],[3,4]];									 inv(a) @ a", { 1,0,0,1 }, { 2,2 }));
+	BP_ASSERT(test_expr("a=[[1,2,3],[4,5,6],[7,8,9]];						 inv(a) @ a", { 1,0,0, 0,1,0, 0,0,1 }, { 3,3 }));
+	BP_ASSERT(test_expr("a=[[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]]; inv(a) @ a", { 1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1 }, { 4,4 }));
+
 	BP_ASSERT(test_expr("abs(-1)+abs(0)+abs(1)", {2}));
 	BP_ASSERT(test_expr("floor(-3.5)", {-4}, {}));
 	BP_ASSERT(test_expr("ceil(-3.5)", {-3}, {}));