Attempt at computing spectral norm using eigenvalues. Requires comparison operators and static cast to double, not possible for Bparser

Author: lukas-valenta-tul

include/array.hh

@@ -1257,11 +1257,17 @@ public:
 		case 2: //matrix
 		{
 			//Spectral norm
-			Throw() << "norm2(matrix) not yet implemented" << "\n";
-			Shape s; //empty Shape for scalar
+			Throw() << "norm2(matrix) is not yet possible" << "\n";
+			/*Shape s; //empty Shape for scalar
 			Array r(s);
-			//r.elements_[0U] = *wrap_array(a);
-			return r;
+
+			Eigen::MatrixX<details::ScalarWrapper> m( wrap_array(a) );
+
+			r.elements_[0U] = *details::sqrt((m.adjoint()*m).eigenvalues().real().maxCoeff());
+			//computing eigenvalues would require static cast to double and comparison operators (<,<=,>,>=,!=,==)
+			//something which we cannot support
+			return r;*/
+			break;
 		}
 		default:
 			Throw() << "Norms are not avaiable for ND tensors" << "\n";

include/scalar_wrapper.hh

@@ -12,6 +12,7 @@
 
 #include "scalar_node.hh"
 #include <Eigen/Core>
+//#include <Eigen/Eigenvalues> //impossible
 
 namespace bparser {
 	namespace details {
@@ -23,6 +24,10 @@ namespace bparser {
 			ScalarWrapper(double d) : node(ScalarNode::create_const(d)) { ; }
 			ScalarWrapper(ScalarNodePtr existing_ptr) : node(existing_ptr) { ; }
 
+			inline ScalarWrapper operator+() const {
+				return ScalarWrapper(*this);
+			}
+
 			inline ScalarWrapper operator-() const {
 				return un_op<_minus_>(*this);
 			}
@@ -36,14 +41,29 @@ namespace bparser {
 				return bin_op<_add_>(*this, b);
 			}
 
+			inline ScalarWrapper& operator-=(const ScalarWrapper& b) {
+				node = bin_op<_sub_>(*this, b).get();
+				return *this;
+			}
+
 			inline ScalarWrapper operator-(const ScalarWrapper& b) const {
 				return bin_op<_sub_>(*this, b);
 			}
 
+			inline ScalarWrapper& operator*=(const ScalarWrapper& b) {
+				node = bin_op<_mul_>(*this, b).get();
+				return *this;
+			}
+
 			inline ScalarWrapper operator*(const ScalarWrapper& b) const {
 				return bin_op<_mul_>(*this, b);
 			}
 
+			inline ScalarWrapper& operator/=(const ScalarWrapper& b) {
+				node = bin_op<_div_>(*this, b).get();
+				return *this;
+			}
+
 			inline ScalarWrapper operator/(const ScalarWrapper& b) const {
 				return bin_op<_div_>(*this, b);
 			}
@@ -53,6 +73,34 @@ namespace bparser {
 					return *(***this).values_ == *(**b).values_;
 				return false;
 			}
+			/* These do not make any sense with what we are trying to achieve
+			inline bool operator!=(const ScalarWrapper& b) const {
+				return !((*this) == b);
+			}
+
+			inline bool operator<(const ScalarWrapper& b) const {
+				if ((*this).is_constant() && (*this).have_same_result_storage(b))
+					return *(***this).values_ < *(**b).values_;
+				return false;
+			}
+
+			inline bool operator<=(const ScalarWrapper& b) const {
+				if ((*this).is_constant() && (*this).have_same_result_storage(b))
+					return *(***this).values_ <= *(**b).values_;
+				return false;
+			}
+
+			inline bool operator>=(const ScalarWrapper& b) const {
+				if ((*this).is_constant() && (*this).have_same_result_storage(b))
+					return *(***this).values_ >= *(**b).values_;
+				return false;
+			}
+
+			inline bool operator>(const ScalarWrapper& b) const {
+				if ((*this).is_constant() && (*this).have_same_result_storage(b))
+					return *(***this).values_ > *(**b).values_;
+				return false;
+			}*/
 
 
 			inline ScalarNodePtr operator*() const { //dereference
@@ -111,7 +159,7 @@ namespace bparser {
 		inline ScalarWrapper abs2(const ScalarWrapper& s) {
 			return s*s;
 		}
-		
+	
 
 		UN_OP(sqrt)
 		//UN_OP(exp)

test/test_parser.cc

@@ -264,7 +264,7 @@ void test_expression() {
 	BP_ASSERT(test_expr("norm1([-4,-3,-2,-1,0,1,2,3,4])", {20}, {}));
 	BP_ASSERT(test_expr("norm1([[-4,-3,-2],[-1,0,1],[2,3,4]])", { 7 }, {}));
 	BP_ASSERT(test_expr("norm2([-4,-3,-2,-1,0,1,2,3,4])", { 7.745966692414834 }, {}));
-	//BP_ASSERT(test_expr("norm2([[-4,-3,-2],[-1,0,1],[2,3,4]])", { 7.3484692283495345 }, {}));
+	//BP_ASSERT(test_expr("norm2([[-4,-3,-2],[-1,0,1],[2,3,4]])", { 7.3484692283495345 }, {})); //Spectral norm uses eigenvalues/singular values. Eigen uses comparison operators in the algorithm. Bparser does not like that
 	BP_ASSERT(test_expr("normfro([[-4,-3,-2],[-1,0,1],[2,3,4]])", { 7.745966692414834 }, {}));
 	BP_ASSERT(test_expr("norminf([-4,-3,-2,-1,0,1,2,3,4])", { 4 }, {}));
 	BP_ASSERT(test_expr("norminf([[-4,-3,-2],[-1,0,1],[2,3,4]])", { 9 }, {}));