// --*-  Mode: C++; c-basic-offset:4; indent-tabs-mode:t; tab-width:4 -*--
// EigenLab
// Version: 1.0.0
// Author: Dr. Marcel Paz Goldschen-Ohm
// Email:  marcel.goldschen@gmail.com
// Copyright (c) 2015 by Dr. Marcel Paz Goldschen-Ohm.
// Licence: MIT
//----------------------------------------

#ifndef EigenLab_H
#define EigenLab_H

#include <Eigen/Dense>
#include <iomanip>
#include <type_traits>
#include <map>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>

// Define both DEBUG and EIGENLAB_DEBUG for step-by-step equation parsing printouts.
#ifndef DEBUG
//#	define DEBUG
#endif

#ifndef EIGENLAB_DEBUG
//#	define EIGENLAB_DEBUG
#endif

#ifdef DEBUG
#	include <iostream>
#endif

namespace EigenLab
{
	//----------------------------------------
	// A wrapper for a matrix whose data is either stored locally or shared.
	//
	// Template typename Derived can be any dynamically sized matrix type supported by Eigen.
	//
	// !!! matrix() promises to ALWAYS return a map to the matrix data whether it's
	// stored locally or externally in some shared memory.
	//
	// !!! local() provides direct access to the local data, but this data is
	// ONLY valid when isLocal() is true. In most cases, you're best off
	// accessing the matrix data via matrix() instead.
	//----------------------------------------
	template <typename Derived = Eigen::MatrixXd>
	class Value
	{
	private:
		// Local matrix data.
		Derived mLocal;
		
		// Map to shared matrix data (map points to mLocal if the data is local).
		// !!! This map promises to ALWAYS point to the matrix data whether it's
		// stored locally in mLocal or externally in some shared memory.
		Eigen::Map<Derived> mShared;
		
		// Flag indicating whether the local data is being used.
		bool mIsLocal;
		
	public:
		// Access mapped data (whether its local or shared).
		// !!! matrix() promises to ALWAYS return a map to the matrix data whether it's
		// stored locally in mLocal or externally in some shared memory.
		inline Eigen::Map<Derived> & matrix() { return mShared; }
		inline const Eigen::Map<Derived> & matrix() const { return mShared; }
		
		// Access local data.
		// !!! WARNING! This data is ONLY valid if isLocal() is true.
		// !!! WARNING! If you change the local data via this method, you MUST call mapLocal() immediately afterwards.
		// In most cases, you're best off accessing the matrix data via matrix() instead.
		inline Derived & local() { return mLocal; }
		inline const Derived & local() const { return mLocal; }
		
		// Is mapped data local?
		inline bool isLocal() const { return mIsLocal; }
		
		// Set mapped data to point to local data.
		inline void mapLocal() { new (& mShared) Eigen::Map<Derived>(mLocal.data(), mLocal.rows(), mLocal.cols()); mIsLocal = true; }
		
		//  Copy shared data to local data (if needed).
		inline void copySharedToLocal() { if(!isLocal()) { mLocal = mShared; mapLocal(); } }
		
		// Set local data.
		Value() : mLocal(1, 1), mShared(mLocal.data(), mLocal.rows(), mLocal.cols()), mIsLocal(true) {}
		Value(const typename Derived::Scalar s) : mLocal(Derived::Constant(1, 1, s)), mShared(mLocal.data(), mLocal.rows(), mLocal.cols()), mIsLocal(true) {}
		Value(const Derived & mat) : mLocal(mat), mShared(mLocal.data(), mLocal.rows(), mLocal.cols()), mIsLocal(true) {}
		inline void setLocal(const typename Derived::Scalar s) { mLocal = Derived::Constant(1, 1, s); mapLocal(); }
		inline void setLocal(const Eigen::MatrixBase<Derived> & mat) { mLocal = mat; mapLocal(); }
		inline void setLocal(const Value & val) { mLocal = val.matrix(); mapLocal(); }
		inline void setLocal(const typename Derived::Scalar * data, size_t rows = 1, size_t cols = 1) { setShared(data, rows, cols); copySharedToLocal(); }
		inline Value & operator = (const typename Derived::Scalar s) { setLocal(s); return (* this); }
		inline Value & operator = (const Derived & mat) { setLocal(mat); return (* this); }
		
		// Set shared data.
		Value(const typename Derived::Scalar * data, size_t rows = 1, size_t cols = 1) : mShared(const_cast<typename Derived::Scalar *>(data), rows, cols), mIsLocal(false) {}
		inline void setShared(const typename Derived::Scalar * data, size_t rows = 1, size_t cols = 1) { new (& mShared) Eigen::Map<Derived>(const_cast<typename Derived::Scalar *>(data), rows, cols); mIsLocal = false; }
		inline void setShared(const Derived & mat) { setShared(mat.data(), mat.rows(), mat.cols()); }
		inline void setShared(const Value & val) { setShared(val.matrix().data(), val.matrix().rows(), val.matrix().cols()); }
		
		// Set to local or shared data dependent on whether val maps its own local data or some other shared data.
		Value(const Value & val) : mLocal(1, 1), mShared(mLocal.data(), mLocal.rows(), mLocal.cols()) { (* this) = val; }
		inline Value & operator = (const Value & val) { if(val.isLocal()) { setLocal(val); } else { setShared(val); } return (* this); }
	};
	typedef Value<Eigen::MatrixXd> ValueXd;
	typedef Value<Eigen::MatrixXf> ValueXf;
	typedef Value<Eigen::MatrixXi> ValueXi;
	
	// check if a class has a comparison operator (ie. std::complex does not)
	template<typename T>
	struct has_operator_lt_impl
	{
		template<class U>
		static auto test(U*) -> decltype(std::declval<U>() < std::declval<U>());
		template<typename>
		static auto test(...) -> std::false_type;
		using type = typename std::is_same<bool, decltype(test<T>(0))>::type;
	};
	template<typename T>
	struct has_operator_lt : has_operator_lt_impl<T>::type {};
	
	//----------------------------------------
	// Equation parser.
	//
	// Template typename Derived can be any dynamically sized matrix type supported by Eigen.
	//----------------------------------------
	template <typename Derived = Eigen::MatrixXd>
	class Parser
	{
	public:
		// A map to hold named values.
		typedef std::map<std::string, Value<Derived> > ValueMap;
		
	private:
		// Variables are stored in a map of named values.
		ValueMap mVariables;
		
		// Operator symbols and function names used by the parser.
		std::string mOperators1, mOperators2;
		std::vector<std::string> mFunctions;
		
		// Expressions are parsed by first splitting them into chunks.
		struct Chunk {
			std::string field;
			int type;
			Value<Derived> value;
			int row0, col0, rows, cols;
			Chunk(const std::string & str = "", int t = -1, const Value<Derived> & val = Value<Derived>()) : field(str), type(t), value(val), row0(-1), col0(-1), rows(-1), cols(-1) {}
		};
		enum ChunkType { VALUE = 0, VARIABLE, OPERATOR, FUNCTION };
		typedef std::vector<Chunk> ChunkArray;
		typedef typename Derived::Index Index;
		bool mCacheChunkedExpressions;
		std::map<std::string, ChunkArray> mCachedChunkedExpressions;
		
	public:
		// Access to named variables.
		// !!! WARNING! var(name) will create the variable name if it does not already exist.
		inline ValueMap & vars() { return mVariables; }
		inline Value<Derived> & var(const std::string & name) { return mVariables[name]; }
		
		// Check if a variable exists.
		inline bool hasVar(const std::string & name) { return isVariable(name); }
		
		// Delete a variable.
		inline void clearVar(const std::string & name) { typename ValueMap::iterator it = mVariables.find(name); if(it != mVariables.end()) mVariables.erase(it); }
		
		// Expression chunk caching.
		inline bool cacheExpressions() const { return mCacheChunkedExpressions; }
		inline void setCacheExpressions(bool b) { mCacheChunkedExpressions = b; }
		inline void clearCachedExpressions() { mCachedChunkedExpressions.clear(); }
		
		Parser();
		~Parser() { clearCachedExpressions(); }
		
		// Evaluate an expression and return the result in a value wrapper.
		Value<Derived> eval(const std::string & expression);
		
	private:
		void splitEquationIntoChunks(const std::string & expression, ChunkArray & chunks, std::string & code);
		std::string::const_iterator findClosingBracket(const std::string & str, const std::string::const_iterator openingBracket, const char closingBracket) const;
		std::vector<std::string> splitArguments(const std::string & str, const char delimeter) const;
		void evalIndexRange(const std::string & str, int * first, int * last, int numIndices);
		void evalMatrixExpression(const std::string & str, Value<Derived> & mat);
		void evalFunction(const std::string & name, std::vector<std::string> & args, Value<Derived> & result);
		bool evalFunction_1_lt(const std::string & name, Value<Derived> & arg, Value<Derived> & result, std::false_type);
		bool evalFunction_1_lt(const std::string & name, Value<Derived> & arg, Value<Derived> & result, std::true_type);
		bool evalFunction_2_lt(const std::string & name, Value<Derived> & arg0, Value<Derived> & arg1, Value<Derived> & result, std::false_type);
		bool evalFunction_2_lt(const std::string & name, Value<Derived> & arg0, Value<Derived> & arg1, Value<Derived> & result, std::true_type);
		
		void evalNumericRange(const std::string & str, Value<Derived> & mat);
		inline bool isVariable(const std::string & name) const { return mVariables.count(name) > 0; }
		inline bool isOperator(const char c) const { return (std::find(mOperators1.begin(), mOperators1.end(), c) != mOperators1.end()); }
		bool isOperator(const std::string & str) const;
		inline bool isFunction(const std::string & str) const { return (std::find(mFunctions.begin(), mFunctions.end(), str) != mFunctions.end()); }
		void evalIndices(ChunkArray & chunks);
		void evalNegations(ChunkArray & chunks);
		void evalPowers(ChunkArray & chunks);
		void evalMultiplication(ChunkArray & chunks);
		void evalAddition(ChunkArray & chunks);
		void evalAssignment(ChunkArray & chunks);
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		void printChunks(ChunkArray & chunks, size_t maxRows = 2, size_t maxCols = 2, int precision = 0);
		void printVars(size_t maxRows = 2, size_t maxCols = 2, int precision = 0);
		std::string textRepresentation(Value<Derived> & val, size_t maxRows = 2, size_t maxCols = 2, int precision = 0);
#	endif
#endif
		
	public:
		static std::string trim(const std::string & str);
		static std::vector<std::string> split(const std::string & str, const char delimeter);
		template <typename T> static bool isNumber(const std::string & str, T * num = 0);
		template <typename T> static T stringToNumber(const std::string & str);
		template <typename T> static std::string numberToString(T num, int precision = 0);
#ifdef DEBUG
		void test_w_lt(size_t & numFails, typename Derived::Scalar & s, Derived & a34, Derived & b34, Derived & c43, Derived & v, std::true_type);
		void test_w_lt(size_t & numFails, typename Derived::Scalar & s, Derived & a34, Derived & b34, Derived & c43, Derived & v, std::false_type);
		size_t test();
#endif
	};
	typedef Parser<Eigen::MatrixXd> ParserXd;
	typedef Parser<Eigen::MatrixXf> ParserXf;
	typedef Parser<Eigen::MatrixXi> ParserXi;
	
	//----------------------------------------
	// Function definitions.
	//----------------------------------------
	template <typename Derived>
	Parser<Derived>::Parser() :
		mOperators1("+-*/^()[]="),
		mOperators2(".+.-.*./.^"),
		mCacheChunkedExpressions(false)
	{
		// Coefficient-wise operations.
		mFunctions.push_back("abs");
		mFunctions.push_back("sqrt");
		mFunctions.push_back("exp");
		mFunctions.push_back("log");
		mFunctions.push_back("log10");
		mFunctions.push_back("sin");
		mFunctions.push_back("cos");
		mFunctions.push_back("tan");
		mFunctions.push_back("asin");
		mFunctions.push_back("acos");

		// Matrix reduction operations.
		mFunctions.push_back("trace");
		mFunctions.push_back("norm");
		mFunctions.push_back("size");
		if (has_operator_lt<typename Derived::Scalar>::value) {
			mFunctions.push_back("min");
			mFunctions.push_back("max");
			mFunctions.push_back("absmax");
		}
		mFunctions.push_back("mean");
		mFunctions.push_back("sum");
		mFunctions.push_back("prod");

		// Matrix operations.
		mFunctions.push_back("transpose");
		mFunctions.push_back("conjugate");
		mFunctions.push_back("adjoint");

		// Matrix initializers.
		mFunctions.push_back("zeros");
		mFunctions.push_back("ones");
		mFunctions.push_back("eye");
	}
	
	template <typename Derived>
	Value<Derived> Parser<Derived>::eval(const std::string & expression)
	{
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		std::cout << "---" << std::endl;
		std::cout << "EXPRESSION: " << expression << std::endl;
#	endif
#endif
		ChunkArray chunks;
		std::string code;
		splitEquationIntoChunks(trim(expression), chunks, code);
		evalIndices(chunks);
		evalNegations(chunks);
		evalPowers(chunks);
		evalMultiplication(chunks);
		evalAddition(chunks);
		evalAssignment(chunks);
		if(chunks.size() != 1)
			throw std::runtime_error("Failed to reduce expression '" + expression + "' to a single value.");
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		std::cout << "---" << std::endl;
#	endif
#endif
		if(chunks[0].type == VARIABLE) {
			if(!isVariable(chunks[0].field))
				throw std::runtime_error("Unknown variable '" + chunks[0].field + "'.");
			return mVariables[chunks[0].field];
		}
		return chunks[0].value;
	}
	
	template <typename Derived>
	void Parser<Derived>::splitEquationIntoChunks(const std::string & expression, ChunkArray & chunks, std::string & code)
	{
		if(cacheExpressions()) {
			if(mCachedChunkedExpressions.count(expression) > 0) {
				chunks = mCachedChunkedExpressions.at(expression);
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
				std::cout << "CACHED CHUNKS: "; printChunks(chunks); std::cout << std::endl;
#	endif
#endif
				return;
			}
		}

		for(std::string::const_iterator it = expression.begin(); it != expression.end();)
		{
			int prevType = (chunks.size() ? chunks.back().type : -1);
			char ci = (* it);
			if(isspace(ci)) {
				// Ignore whitespace.
				it++;
			} else if(ci == '(' && (prevType == VALUE || prevType == VARIABLE)) {
				// Index group.
				std::string::const_iterator jt = findClosingBracket(expression, it, ')');
				if(jt == expression.end())
					throw std::runtime_error("Missing closing bracket for '" + std::string(it, jt) + "'.");
				std::string field = std::string(it + 1, jt); // Outer brackets stripped.
				if(prevType == VARIABLE) {
					if(!isVariable(chunks.back().field))
						throw std::runtime_error("Unknown variable '" + chunks.back().field + "'.");
					chunks.back().value.setShared(var(chunks.back().field));
				}
				int first, last;
				int rows = int(chunks.back().value.matrix().rows());
				int cols = int(chunks.back().value.matrix().cols());
				std::vector<std::string> args = splitArguments(field, ',');
				if(args.size() == 1) {
					if(cols == 1) {
						evalIndexRange(args[0], & first, & last, rows);
						chunks.back().row0 = first;
						chunks.back().col0 = 0;
						chunks.back().rows = last + 1 - first;//(last == -1 ? int(chunks.back().value.matrix().rows()) : last + 1) - first;
						chunks.back().cols = 1;
					} else if(rows == 1) {
						evalIndexRange(args[0], & first, & last, cols);
						chunks.back().row0 = 0;
						chunks.back().col0 = first;
						chunks.back().rows = 1;
						chunks.back().cols = last + 1 - first;//(last == -1 ? int(chunks.back().value.matrix().cols()) : last + 1) - first;
					} else {
						throw std::runtime_error("Missing row or column indices for '(" + chunks.back().field + "(" + field + ")'.");
					}
				} else if(args.size() == 2) {
					evalIndexRange(args[0], & first, & last, rows);
					chunks.back().row0 = first;
					chunks.back().rows = last + 1 - first;//(last == -1 ? int(chunks.back().value.matrix().rows()) : last + 1) - first;
					evalIndexRange(args[1], & first, & last, cols);
					chunks.back().col0 = first;
					chunks.back().cols = last + 1 - first;//(last == -1 ? int(chunks.back().value.matrix().cols()) : last + 1) - first;
				} else {
					throw std::runtime_error("Invalid index expression '" + chunks.back().field + "(" + field + ")'.");
				}
				code += "i";
				it = jt + 1;
			} else if(ci == '(' && prevType == FUNCTION) {
				// Function argument group.
				std::string::const_iterator jt = findClosingBracket(expression, it, ')');
				if(jt == expression.end())
					throw std::runtime_error("Missing closing bracket for '" + std::string(it, jt) + "'.");
				std::string field = std::string(it + 1, jt); // Outer brackets stripped.
				std::vector<std::string> args = splitArguments(field, ',');
				evalFunction(chunks.back().field, args, chunks.back().value);
				chunks.back().field += "(" + field + ")";
				chunks.back().type = VALUE;
				code += (chunks.back().value.matrix().size() == 1 ? "n" : "M");
				it = jt + 1;
			} else if(ci == '(') {
				// Recursively evaluate group to a single value.
				std::string::const_iterator jt = findClosingBracket(expression, it, ')');
				if(jt == expression.end())
					throw std::runtime_error("Missing closing bracket for '" + std::string(it, jt) + "'.");
				std::string field = std::string(it + 1, jt); // Outer brackets stripped.
				chunks.push_back(Chunk(field, prevType = VALUE, eval(field)));
				code += (chunks.back().value.matrix().size() == 1 ? "n" : "M");
				it = jt + 1;
			} else if(ci == '[') {
				// Evaluate matrix.
				if(prevType == VALUE || prevType == VARIABLE)
					throw std::runtime_error("Invalid operation '" + chunks.back().field + std::string(1, ci) + "'.");
				std::string::const_iterator jt = findClosingBracket(expression, it, ']');
				if(jt == expression.end())
					throw std::runtime_error("Missing closing bracket for '" + std::string(it, jt) + "'.");
				std::string field = std::string(it + 1, jt); // Outer brackets stripped.
				chunks.push_back(Chunk("[" + field + "]", prevType = VALUE));
				evalMatrixExpression(field, chunks.back().value);
				code += (chunks.back().value.matrix().size() == 1 ? "n" : "M");
				it = jt + 1;
			} else if(it + 1 != expression.end() && isOperator(std::string(it, it + 2))) {
				// Double character operator.
				std::string field = std::string(it, it + 2);
				chunks.push_back(Chunk(field, prevType = OPERATOR));
				code += field;
				it += 2;
			} else if(isOperator(ci)) {
				// Single character operator.
				std::string field = std::string(1, ci);
				chunks.push_back(Chunk(field, prevType = OPERATOR));
				code += field;
				it++;
			} else {
				// Non-operator: value range, number, function, or variable name.
				std::string::const_iterator jt = it + 1;
				// accept fp-strings, ie [+-]
				unsigned char state = 1;
				for(std::string::const_iterator kt = it; state && kt != expression.end(); kt++) {
					unsigned char token;
					if (*kt == ' ') token = 0;
					else if (*kt == '+' || *kt == '-') token = 1;
					else if (isdigit(*kt)) token = 2;
					else if (*kt == '.') token = 3;
					else if (*kt == 'e' || *kt == 'E') token = 4;
					else break;
					const static char nextstate[9][5] = {{0},
														 {1, 2, 3, 4, 0},
														 {0, 0, 3, 4, 0},
														 {0, 0, 3, 5, 6},
														 {0, 0, 5, 0, 0},
														 {0, 0, 5, 0, 6},
														 {0, 7, 8, 0, 0},
														 {0, 0, 8, 0, 0},
														 {0, 0, 8, 0, 0}};
					//WARN("state=" << (int)state << " token(" << *kt << ")=" << (int)token
					//	 << " nextstate = " << (int)nextstate[state][token] << "\n");
					state = nextstate[state][token];
					if (state == 8) jt = kt;
				}
				for(; jt != expression.end(); jt++) {
					if(isOperator(*jt) || (jt + 1 != expression.end() && isOperator(std::string(jt, jt + 2))))
						break;
				}
				std::string field = trim(std::string(it, jt));
				if(prevType == VALUE || prevType == VARIABLE)
					throw std::runtime_error("Invalid operation '" + chunks.back().field + field + "'.");
				double num;
				if(field.find(":") != std::string::npos) {
					// Numeric range.
					chunks.push_back(Chunk(field, prevType = VALUE));
					evalNumericRange(field, chunks.back().value);
					code += (chunks.back().value.matrix().size() == 1 ? "n" : "M");
				} else if(isNumber<double>(field, & num)) {
					// Number.
					chunks.push_back(Chunk(field, prevType = VALUE, Value<Derived>(num)));
					code += "n";
				} else if(isVariable(field)) {
					// Local variable.
					chunks.push_back(Chunk(field, prevType = VARIABLE));
					code += (mVariables[field].matrix().size() == 1 ? "vn" : "vM");
				} else if(isFunction(field)) {
					// Function.
					chunks.push_back(Chunk(field, prevType = FUNCTION));
				} else {
					// New undefined variable.
					chunks.push_back(Chunk(field, prevType = VARIABLE));
					code += "vU";
				}
				it = jt;
			}
		} // it
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		std::cout << "CHUNKS: "; printChunks(chunks); std::cout << std::endl;
		std::cout << "CODE: " << code << std::endl;
#	endif
#endif
		if(cacheExpressions())
			mCachedChunkedExpressions[expression] = chunks;
	}
	
	template <typename Derived>
	std::string::const_iterator Parser<Derived>::findClosingBracket(const std::string & str, const std::string::const_iterator openingBracket, const char closingBracket) const
	{
		int depth = 1;
		std::string::const_iterator it = openingBracket + 1;
		for(; it != str.end(); it++) {
			if((*it) == (*openingBracket)) depth++;
			else if((*it) == closingBracket) depth--;
			if(depth == 0) return it;
		}
		return str.end();
	}

	template <typename Derived>
	std::vector<std::string> Parser<Derived>::splitArguments(const std::string & str, const char delimeter) const
	{
		std::vector<std::string> args;
		std::string::const_iterator i0 = str.begin();
		for(std::string::const_iterator it = str.begin(); it != str.end(); it++) {
			if((*it) == '(') it = findClosingBracket(str, it, ')');
			else if((*it) == '[') it = findClosingBracket(str, it, ']');
			else if((*it) == delimeter) {
				args.push_back(trim(std::string(i0, it)));
				i0 = it + 1;
			}
		}
		args.push_back(std::string(i0, str.end()));
		return args;
	}

	template <typename Derived>
	void Parser<Derived>::evalIndexRange(const std::string & str, int * first, int * last, int numIndices)
	{
		if(str.empty())
			throw std::runtime_error("Empty index range.");
		ValueXi valuei;
		ParserXi parseri;
		size_t pos;
		for(std::string::const_iterator it = str.begin(); it != str.end(); it++) {
			if((*it) == ':') {
				std::string firstStr = trim(std::string(str.begin(), it));
				std::string lastStr = trim(std::string(it + 1, str.end()));
				if(firstStr.empty() && lastStr.empty()) {
					(* first) = 0;
					(* last) = numIndices - 1;
					return;
				}
				if(firstStr.empty() || lastStr.empty())
					throw std::runtime_error("Missing indices for '" + str + "'.");
				
				pos = firstStr.find("end");
				if(pos != std::string::npos) {
					firstStr = firstStr.substr(0, pos) + numberToString<int>(numIndices - 1) + firstStr.substr(pos + 3);
				}
				pos = lastStr.find("end");
				if(pos != std::string::npos) {
					lastStr = lastStr.substr(0, pos) + numberToString<int>(numIndices - 1) + lastStr.substr(pos + 3);
				}
				
				valuei = parseri.eval(firstStr);
				if(valuei.matrix().size() != 1)
					throw std::runtime_error("Invalid indices '" + str + "'.");
				(* first) = valuei.matrix()(0, 0);
				
				valuei = parseri.eval(lastStr);
				if(valuei.matrix().size() != 1)
					throw std::runtime_error("Invalid indices '" + str + "'.");
				(* last) = valuei.matrix()(0, 0);
				
				return;
			}
		}
		std::string firstStr = str;
		
		pos = firstStr.find("end");
		if(pos != std::string::npos) {
			firstStr = firstStr.substr(0, pos) + numberToString<int>(numIndices - 1) + firstStr.substr(pos + 3);
		}
		
		valuei = parseri.eval(firstStr);
		if(valuei.matrix().size() != 1)
			throw std::runtime_error("Invalid index '" + str + "'.");
		(* first) = valuei.matrix()(0, 0);
		(* last) = (* first);
	}

	template <typename Derived>
	void Parser<Derived>::evalMatrixExpression(const std::string & str, Value<Derived> & mat)
	{
		// !!! Expression may NOT include outer brackets, although brackets for individual rows are OK.
		std::vector<std::string> rows = splitArguments(str, ';');
		std::vector<std::vector<typename Derived::Scalar> > temp;
		Value<Derived> submatrix;
		size_t row0 = 0, col0 = 0, nrows = 0, ncols = 0;
		size_t pos;
		for(size_t i = 0; i < rows.size(); i++) {
			// Strip row brackets if they exist.
			if(rows[i][0] == '[' && rows[i].back() == ']') rows[i] = rows[i].substr(1, int(rows[i].size()) - 2);
			std::vector<std::string> cols = splitArguments(rows[i], ',');
			col0 = 0;
			ncols = 0;
			for(size_t j = 0; j < cols.size(); j++) {
				pos = cols[j].find(":");
				if(pos != std::string::npos) {
					std::string firstStr = cols[j].substr(0, pos);
					std::string lastStr = cols[j].substr(pos + 1);
					pos = lastStr.find(":");
					if(pos != std::string::npos) {
						std::string stepStr = lastStr.substr(0, pos);
						lastStr = lastStr.substr(pos + 1);
						if(lastStr.find(":") != std::string::npos)
							throw std::runtime_error("Invalid matrix definition '[" + str + "]'. Invalid range '" + cols[j] + "'.");
						// first:step:last
						Value<Derived> first = eval(firstStr);
						Value<Derived> step = eval(stepStr);
						Value<Derived> last = eval(lastStr);
						if(first.matrix().size() != 1 || step.matrix().size() != 1 || last.matrix().size() != 1)
							throw std::runtime_error("Invalid matrix definition '[" + str + "]'. Invalid range '" + cols[j] + "'.");
						typename Derived::RealScalar sfirst = std::real(first.matrix()(0));
						typename Derived::RealScalar sstep = std::real(step.matrix()(0));
						typename Derived::RealScalar slast = std::real(last.matrix()(0));
						if(sfirst == slast) {
							submatrix.local().setConstant(1, 1, sfirst);
							submatrix.mapLocal();
						} else if((slast - sfirst >= 0 && sstep > 0) || (slast - sfirst <= 0 && sstep < 0)) {
							int n = floor((slast - sfirst) / sstep) + 1;
							submatrix.local().resize(1, n);
							for(int k = 0; k < n; ++k)
								submatrix.local()(0, k) = sfirst + k * sstep;
							submatrix.mapLocal();
						} else {
							throw std::runtime_error("Invalid matrix definition '[" + str + "]'. Invalid range '" + cols[j] + "'.");
						}
					} else {
						// first:last => first:1:last
						Value<Derived> first = eval(firstStr);
						Value<Derived> last = eval(lastStr);
						if(first.matrix().size() != 1 || last.matrix().size() != 1)
							throw std::runtime_error("Invalid matrix definition '[" + str + "]'. Invalid range '" + cols[j] + "'.");
						typename Derived::RealScalar sfirst = std::real(first.matrix()(0));
						typename Derived::RealScalar slast = std::real(last.matrix()(0));
						if(sfirst == slast) {
							submatrix.local().setConstant(1, 1, sfirst);
							submatrix.mapLocal();
						} else if(slast - sfirst >= 0) {
							int n = floor(slast - sfirst) + 1;
							submatrix.local().resize(1, n);
							for(int k = 0; k < n; ++k)
								submatrix.local()(0, k) = sfirst + k;
							submatrix.mapLocal();
						} else {
							throw std::runtime_error("Invalid matrix definition '[" + str + "]'. Invalid range '" + cols[j] + "'.");
						}
					}
				} else {
					submatrix = eval(cols[j]);
				}
				if(j > 0 && size_t(submatrix.matrix().cols()) != nrows)
					throw std::runtime_error("Invalid matrix definition '[" + str + "]'. Successive column entries '" + cols[int(j) - 1] + "' and '" + cols[j] + "' do not have the same number of rows.");
				nrows = submatrix.matrix().rows();
				ncols += submatrix.matrix().cols();
				temp.resize(row0 + submatrix.matrix().rows());
				for(size_t row = 0; row < size_t(submatrix.matrix().rows()); row++) {
					temp[row0 + row].resize(col0 + submatrix.matrix().cols());
					for(size_t col = 0; col < size_t(submatrix.matrix().cols()); col++)
						temp[row0 + row][col0 + col] = submatrix.matrix()(row, col);
				}
				col0 += submatrix.matrix().cols();
			}
			if(row0 > 0 && ncols != temp[int(row0) - 1].size())
				throw std::runtime_error("Invalid matrix definition '[" + str + "]'. Successive row entries '" + rows[int(i) - 1] + "' and '" + rows[i] + "' do not have the same number of columns.");
			row0 += nrows;
		}
		nrows = temp.size();
		if(nrows == 0) return;
		ncols = temp[0].size();
		mat.setLocal(Derived(nrows, ncols));
		for(size_t row = 0; row < nrows; row++) {
			for(size_t col = 0; col < ncols; col++)
				mat.matrix()(row, col) = temp[row][col];
		}
		mat.mapLocal();
	}

	template <typename Derived>
	bool Parser<Derived>::evalFunction_1_lt(const std::string & name, Value<Derived> & arg, Value<Derived> & result, std::true_type)
	{
		if(name == "min") {
			result.setLocal(arg.matrix().minCoeff());
			return true;
		} else if(name == "max") {
			result.setLocal(arg.matrix().maxCoeff());
			return true;
		} else if(name == "absmax") {
			typename Derived::Scalar minimum = arg.matrix().minCoeff();
			typename Derived::Scalar maximum = arg.matrix().maxCoeff();
			result.setLocal(std::abs(maximum) >= std::abs(minimum) ? maximum : minimum);
			return true;
		}
		return false;
	}

	template <typename Derived>
	bool Parser<Derived>::evalFunction_1_lt(const std::string &/*name*/, Value<Derived> &/*arg0*/, Value<Derived> &/*result*/, std::false_type)
	{
		return false;
	}

	template <typename Derived>
	bool Parser<Derived>::evalFunction_2_lt(const std::string & name, Value<Derived> & arg0, Value<Derived> & arg1, Value<Derived> & result, std::true_type)
	{
        if(name == "min") {
            if(arg1.matrix().size() != 1)
                throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
            int dim = floor(std::real(arg1.matrix()(0, 0)));
            if((dim != 0 && dim != 1) || dim != std::real(arg1.matrix()(0, 0)))
                throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
			if(dim == 0) {
				result.local() = arg0.matrix().colwise().minCoeff();
				result.mapLocal();
				return true;
			} else if(dim == 1) {
				result.local() = arg0.matrix().rowwise().minCoeff();
				result.mapLocal();
				return true;
			}
        } else if(name == "max") {
            if(arg1.matrix().size() != 1)
                throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
            int dim = floor(std::real(arg1.matrix()(0, 0)));
            if((dim != 0 && dim != 1) || dim != std::real(arg1.matrix()(0, 0)))
                throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
			if(dim == 0) {
				result.local() = arg0.matrix().colwise().maxCoeff();
				result.mapLocal();
				return true;
			} else if(dim == 1) {
				result.local() = arg0.matrix().rowwise().maxCoeff();
				result.mapLocal();
				return true;
			}
        } else if(name == "absmax") {
            if(arg1.matrix().size() != 1)
                throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
            int dim = floor(std::real(arg1.matrix()(0, 0)));
            if((dim != 0 && dim != 1) || dim != std::real(arg1.matrix()(0, 0)))
                throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
			if(dim == 0) {
				result.local() = arg0.matrix().colwise().maxCoeff();
				result.mapLocal();
				Derived minimum = arg0.matrix().colwise().minCoeff();
				for(size_t i = 0; i < size_t(result.matrix().size()); i++) {
					if(std::abs(result.matrix()(i)) < std::abs(minimum(i)))
						result.matrix()(i) = minimum(i);
				}
				return true;
			} else if(dim == 1) {
				result.local() = arg0.matrix().rowwise().maxCoeff();
				result.mapLocal();
				Derived minimum = arg0.matrix().rowwise().minCoeff();
				for(size_t i = 0; i < size_t(result.matrix().size()); i++) {
					if(std::abs(result.matrix()(i)) < std::abs(minimum(i)))
						result.matrix()(i) = minimum(i);
				}
				return true;
            }
		}
		return false;
	}

	template <typename Derived>
	bool Parser<Derived>::evalFunction_2_lt(const std::string &/*name*/, Value<Derived> &/*arg0*/, Value<Derived> &/*arg1*/, Value<Derived> &/*result*/, std::false_type)
	{
		return false;
	}
	
	template <typename Derived>
	void Parser<Derived>::evalFunction(const std::string & name, std::vector<std::string> & args, Value<Derived> & result)
	{
		if(args.size() == 1) {
			Value<Derived> arg = eval(args[0]);
			if(name == "abs") {
				result.local() = arg.matrix().array().abs().template cast<typename Derived::Scalar>();
				result.mapLocal();
				return;
			} else if(name == "sqrt") {
				result.local() = arg.matrix().array().sqrt();
				result.mapLocal();
				return;
			} else if(name == "exp") {
				result.local() = arg.matrix().array().exp();
				result.mapLocal();
				return;
			} else if(name == "log") {
				result.local() = arg.matrix().array().log();
				result.mapLocal();
				return;
			} else if(name == "log10") {
                result.local() = arg.matrix().array().log();
                result.local() *= (1.0 / log(10));
				result.mapLocal();
				return;
			} else if(name == "sin") {
				result.local() = arg.matrix().array().sin();
				result.mapLocal();
				return;
			} else if(name == "cos") {
				result.local() = arg.matrix().array().cos();
				result.mapLocal();
				return;
			} else if(name == "tan") {
				result.local() = arg.matrix().array().tan();
				result.mapLocal();
				return;
			} else if(name == "acos") {
				result.local() = arg.matrix().array().acos();
				result.mapLocal();
				return;
			} else if(name == "asin") {
				result.local() = arg.matrix().array().asin();
				result.mapLocal();
				return;
			} else if(name == "trace") {
				result.setLocal(arg.matrix().trace());
				return;
			} else if(name == "norm") {
				result.setLocal(arg.matrix().norm());
				return;
			} else if (evalFunction_1_lt(name, arg, result, has_operator_lt<typename Derived::Scalar>())) {
				return;
			} else if(name == "mean") {
				result.setLocal(arg.matrix().mean());
				return;
			} else if(name == "sum") {
				result.setLocal(arg.matrix().sum());
				return;
			} else if(name == "prod") {
				result.setLocal(arg.matrix().prod());
				return;
			} else if(name == "transpose") {
				result.local() = arg.matrix().transpose();
				result.mapLocal();
				return;
			} else if(name == "conjugate") {
				result.local() = arg.matrix().conjugate();
				result.mapLocal();
				return;
			} else if(name == "adjoint") {
				result.local() = arg.matrix().adjoint();
				result.mapLocal();
				return;
            } else if(name == "zeros") {
                if(arg.matrix().size() != 1)
                    throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + ")'.");
                int N = floor(std::real(arg.matrix()(0, 0)));
                if(N <= 0 || N != std::real(arg.matrix()(0, 0)))
                    throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + ")'.");
                result.local() = Derived::Zero(N, N);
                result.mapLocal();
                return;
            } else if(name == "ones") {
                if(arg.matrix().size() != 1)
                    throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + ")'.");
                int N = floor(std::real(arg.matrix()(0, 0)));
                if(N <= 0 || N != std::real(arg.matrix()(0, 0)))
                    throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + ")'.");
                result.local() = Derived::Ones(N, N);
                result.mapLocal();
                return;
            } else if(name == "eye") {
                if(arg.matrix().size() != 1)
                    throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + ")'.");
                int N = floor(std::real(arg.matrix()(0, 0)));
                if(N <= 0 || N != std::real(arg.matrix()(0, 0)))
                    throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + ")'.");
                result.local() = Derived::Identity(N, N);
                result.mapLocal();
                return;
			}
			else {
				throw std::runtime_error("Invalid function '" + name + "(" + args[0] + ")'.");
			}
		} else if(args.size() == 2) {
			Value<Derived> arg0 = eval(args[0]);
			Value<Derived> arg1 = eval(args[1]);
            if(name == "size") {
                if(arg1.matrix().size() != 1)
                    throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
                int dim = floor(std::real(arg1.matrix()(0, 0)));
                if((dim != 0 && dim != 1) || dim != std::real(arg1.matrix()(0, 0)))
                    throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
				if(dim == 0) {
					result.setLocal((typename Derived::Scalar) arg0.matrix().rows());
					return;
				} else if(dim == 1) {
					result.setLocal((typename Derived::Scalar) arg0.matrix().cols());
					return;
                }
			} else if (evalFunction_2_lt(name, arg0, arg1, result, has_operator_lt<typename Derived::Scalar>())) {
				return;
			} else if(name == "mean") {
                if(arg1.matrix().size() != 1)
                    throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
                int dim = floor(std::real(arg1.matrix()(0, 0)));
                if((dim != 0 && dim != 1) || dim != std::real(arg1.matrix()(0, 0)))
                    throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
				if(dim == 0) {
					result.local() = arg0.matrix().colwise().mean();
					result.mapLocal();
					return;
				} else if(dim == 1) {
					result.local() = arg0.matrix().rowwise().mean();
					result.mapLocal();
					return;
				}
            } else if(name == "sum") {
                if(arg1.matrix().size() != 1)
                    throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
                int dim = floor(std::real(arg1.matrix()(0, 0)));
                if((dim != 0 && dim != 1) || dim != std::real(arg1.matrix()(0, 0)))
                    throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
				if(dim == 0) {
					result.local() = arg0.matrix().colwise().sum();
					result.mapLocal();
					return;
				} else if(dim == 1) {
					result.local() = arg0.matrix().rowwise().sum();
					result.mapLocal();
					return;
				}
            } else if(name == "prod") {
                if(arg1.matrix().size() != 1)
                    throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
                int dim = floor(std::real(arg1.matrix()(0, 0)));
                if((dim != 0 && dim != 1) || dim != std::real(arg1.matrix()(0, 0)))
                    throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
				if(dim == 0) {
					result.local() = arg0.matrix().colwise().prod();
					result.mapLocal();
					return;
				} else if(dim == 1) {
					result.local() = arg0.matrix().rowwise().prod();
					result.mapLocal();
					return;
                }
            } else if(name == "zeros") {
                if((arg0.matrix().size() != 1) || (arg1.matrix().size() != 1))
                    throw std::runtime_error("Invalid dimension arguments for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
                int rows = floor(std::real(arg0.matrix()(0, 0)));
                int cols = floor(std::real(arg1.matrix()(0, 0)));
                if(rows <= 0 || cols <= 0 || rows != std::real(arg0.matrix()(0, 0)) || cols != std::real(arg1.matrix()(0, 0)))
                    throw std::runtime_error("Invalid dimension arguments for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
                result.local() = Derived::Zero(rows, cols);
                result.mapLocal();
                return;
            } else if(name == "ones") {
                if((arg0.matrix().size() != 1) || (arg1.matrix().size() != 1))
                    throw std::runtime_error("Invalid dimension arguments for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
                int rows = floor(std::real(arg0.matrix()(0, 0)));
                int cols = floor(std::real(arg1.matrix()(0, 0)));
                if(rows <= 0 || cols <= 0 || rows != std::real(arg0.matrix()(0, 0)) || cols != std::real(arg1.matrix()(0, 0)))
                    throw std::runtime_error("Invalid dimension arguments for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
                result.local() = Derived::Ones(rows, cols);
                result.mapLocal();
                return;
            } else if(name == "eye") {
                if((arg0.matrix().size() != 1) || (arg1.matrix().size() != 1))
                    throw std::runtime_error("Invalid dimension arguments for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
                int rows = floor(std::real(arg0.matrix()(0, 0)));
                int cols = floor(std::real(arg1.matrix()(0, 0)));
                if(rows <= 0 || cols <= 0 || rows != std::real(arg0.matrix()(0, 0)) || cols != std::real(arg1.matrix()(0, 0)))
                    throw std::runtime_error("Invalid dimension arguments for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
                result.local() = Derived::Identity(rows, cols);
                result.mapLocal();
                return;
			} else {
				throw std::runtime_error("Invalid function '" + name + "(" + args[0] + "," + args[1] + ")'.");
			}
		}
		std::string argsStr = "(";
		for(size_t i = 0; i < args.size(); i++) {
			if(i > 0) argsStr += ",";
			argsStr += args[i];
		}
		argsStr += ")";
		throw std::runtime_error("Invalid function/arguments for '" + name + argsStr + "'.");
	}
	
	template <typename Derived>
	void Parser<Derived>::evalNumericRange(const std::string & str, Value<Derived> & mat)
	{
		size_t pos = str.find(":");
		if(pos == std::string::npos)
			throw std::runtime_error("Invalid numeric range '" + str + "'.");
		size_t pos2 = str.substr(pos + 1).find(":");
		if(pos2 == std::string::npos) {
			// first:last
			std::string firstStr = str.substr(0, pos);
			std::string lastStr = str.substr(pos + 1);
			Value<Derived> first = eval(firstStr);
			Value<Derived> last = eval(lastStr);
			if(first.matrix().size() != 1 || last.matrix().size() != 1)
				throw std::runtime_error("Invalid numeric range '" + str + "'.");
			typename Derived::RealScalar sfirst = std::real(first.matrix()(0,0));
			typename Derived::RealScalar slast = std::real(last.matrix()(0,0));
			if(sfirst > slast)
				throw std::runtime_error("Invalid numeric range '" + str + "'. Must not reverse.");
			int n = 1 + floor(slast - sfirst);
			mat.local().resize(1, n);
			for(int i = 0; i < n; i++)
				mat.local()(0, i) = sfirst + i;
			mat.mapLocal();
		} else {
			// first:step:last
			pos2 += pos + 1;
			std::string firstStr = str.substr(0, pos);
			std::string stepStr = str.substr(pos + 1, pos2 - pos - 1);
			std::string lastStr = str.substr(pos2 + 1);
			Value<Derived> first = eval(firstStr);
			Value<Derived> step = eval(stepStr);
			Value<Derived> last = eval(lastStr);
			if(first.matrix().size() != 1 || step.matrix().size() != 1 || last.matrix().size() != 1)
				throw std::runtime_error("Invalid numeric range '" + str + "'.");
			typename Derived::RealScalar sfirst = std::real(first.matrix()(0, 0));
			typename Derived::RealScalar sstep = std::real(step.matrix()(0, 0));
			typename Derived::RealScalar slast = std::real(last.matrix()(0, 0));
			if(sfirst == slast) {
				mat = sfirst;
			} else if(sfirst < slast && sstep > 0) {
				int n = 1 + floor((slast - sfirst) / sstep);
				mat.local().resize(1, n);
				for(int i = 0; i < n; i++)
					mat.local()(0, i) = sfirst + i * sstep;
				mat.mapLocal();
			} else if(sfirst > slast && sstep < 0) {
				int n = 1 + floor((slast - sfirst) / sstep);
				mat.local().resize(1, n);
				for(int i = 0; i < n; i++)
					mat.local()(0, i) = sfirst + i * sstep;
				mat.mapLocal();
			} else {
				throw std::runtime_error("Invalid numeric range '" + str + "'.");
			}
		}
	}
	
	template <typename Derived>
	bool Parser<Derived>::isOperator(const std::string & str) const
	{
		if(str.size() == 1)
			return isOperator(str[0]);
		else if(str.size() == 2) {
			size_t pos = mOperators2.find(str);
			return (pos != std::string::npos && pos % 2 == 0);
		}
		return false;
	}
	
	template <typename Derived>
	void Parser<Derived>::evalIndices(ChunkArray & chunks)
	{
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		bool operationPerformed = false;
#	endif
#endif
		for(typename ChunkArray::iterator it = chunks.begin(); it != chunks.end(); it++) {
			if(it->row0 != -1 && (it->type == VALUE || (it->type == VARIABLE && (it + 1 == chunks.end() || (it + 1)->type != OPERATOR || (it + 1)->field != "=")))) {
				if(it->type == VALUE) {
					Derived temp = it->value.local().block(it->row0, it->col0, it->rows, it->cols);
					it->value.local() = temp;
					it->value.mapLocal();
				} else { //if(it->type == VARIABLE) {
					if(!isVariable(it->field))
						throw std::runtime_error("Attempted indexing into uninitialized variable '" + it->field + "'.");
					it->value.local() = mVariables[it->field].matrix().block(it->row0, it->col0, it->rows, it->cols);
					it->value.mapLocal();
					it->type = VALUE;
				}
				it->row0 = -1;
				it->col0 = -1;
				it->rows = -1;
				it->cols = -1;
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
				operationPerformed = true;
#	endif
#endif
			}
		}
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		if(operationPerformed) { std::cout << "i: "; printChunks(chunks); std::cout << std::endl; }
#	endif
#endif
	}
	
	template <typename Derived>
	void Parser<Derived>::evalNegations(ChunkArray & chunks)
	{
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		bool operationPerformed = false;
#	endif
#endif
		if(chunks.size() < 2) return;
		typename ChunkArray::iterator lhs = chunks.begin(), op = chunks.begin(), rhs = op + 1;
		for(; lhs != chunks.end() && op != chunks.end() && rhs != chunks.end();) {
			if(op->type == OPERATOR && op->field == "-" && (op == chunks.begin() || (lhs->type != VALUE && lhs->type != VARIABLE)) && (rhs->type == VALUE || rhs->type == VARIABLE)) {
				if(rhs->type == VALUE)
					rhs->value.matrix().array() *= -1;
				else if(rhs->type == VARIABLE) {
					if(!isVariable(rhs->field))
						throw std::runtime_error("Attempted operation '" + op->field + rhs->field + "' on uninitialized variable '" + rhs->field + "'.");
					rhs->value.local() = mVariables[rhs->field].matrix().array() * -1;
					rhs->value.mapLocal();
					rhs->type = VALUE;
				}
				lhs = chunks.erase(op);
				op = (lhs != chunks.end()) ? lhs + 1 : lhs;
				rhs = (op != chunks.end()) ? op + 1 : op;
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
				operationPerformed = true;
#	endif
#endif
			} else {
				lhs = op;
				op = rhs;
				rhs++;
			}
		}
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		if(operationPerformed) { std::cout << "-: "; printChunks(chunks); std::cout << std::endl; }
#	endif
#endif
	}

	template <typename Derived>
	void Parser<Derived>::evalPowers(ChunkArray & chunks)
	{
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		bool operationPerformed = false;
#	endif
#endif
		if(chunks.size() < 3) return;
		typename ChunkArray::iterator lhs = chunks.begin(), op = lhs + 1, rhs = op + 1;
		for(; lhs != chunks.end() && op != chunks.end() && rhs != chunks.end();) {
			if((op->type == OPERATOR) && (op->field == "^" || op->field == ".^")) {
				if(lhs->type == VARIABLE) {
					if(!isVariable(lhs->field))
						throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + lhs->field + "'.");
					lhs->value.setShared(mVariables[lhs->field]);
				}
				if(rhs->type == VARIABLE) {
					if(!isVariable(rhs->field))
						throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + rhs->field + "'.");
					rhs->value.setShared(mVariables[rhs->field]);
				}
				if(rhs->value.matrix().size() == 1) {
					lhs->value.local() = lhs->value.matrix().array().pow(rhs->value.matrix()(0, 0));
					lhs->value.mapLocal();
					lhs->type = VALUE;
				} else if(lhs->value.matrix().size() == 1) {
					typename Derived::Scalar temp = lhs->value.matrix()(0, 0);
					lhs->value.local().resize(rhs->value.matrix().rows(), rhs->value.matrix().cols());
					for(size_t row = 0; row < size_t(rhs->value.matrix().rows()); row++) {
						for(size_t col = 0; col < size_t(rhs->value.matrix().cols()); col++)
							lhs->value.local()(row, col) = pow(temp, rhs->value.matrix()(row, col));
					}
					lhs->value.mapLocal();
					lhs->type = VALUE;
				} else if(op->field == ".^" && lhs->value.matrix().rows() == rhs->value.matrix().rows() && lhs->value.matrix().cols() == rhs->value.matrix().cols()) {
					lhs->value.local().resize(rhs->value.matrix().rows(), rhs->value.matrix().cols());
					for(size_t row = 0; row < size_t(rhs->value.matrix().rows()); row++) {
						for(size_t col = 0; col < size_t(rhs->value.matrix().cols()); col++)
							lhs->value.local()(row, col) = pow(lhs->value.matrix()(row, col), rhs->value.matrix()(row, col));
					}
					lhs->value.mapLocal();
					lhs->type = VALUE;
				} else {
					throw std::runtime_error("Invalid operand dimensions for operation '" + lhs->field + op->field + rhs->field + "'.");
				}
				chunks.erase(op, rhs + 1);
				op = lhs + 1;
				rhs = (op != chunks.end()) ? op + 1 : op;
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
				operationPerformed = true;
#	endif
#endif
			} else {
				lhs = op;
				op = rhs;
				rhs++;
			}
		}
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		if(operationPerformed) { std::cout << "^: "; printChunks(chunks); std::cout << std::endl; }
#	endif
#endif
	}

	template <typename Derived>
	void Parser<Derived>::evalMultiplication(ChunkArray & chunks)
	{
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		bool operationPerformed = false;
#	endif
#endif
		if(chunks.size() < 3) return;
		typename ChunkArray::iterator lhs = chunks.begin(), op = lhs + 1, rhs = op + 1;
		for(; lhs != chunks.end() && op != chunks.end() && rhs != chunks.end();) {
			if((op->type == OPERATOR) && (op->field == "*" || op->field == "/" || op->field == ".*" || op->field == "./")) {
				if(lhs->type == VARIABLE) {
					if(!isVariable(lhs->field))
						throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + lhs->field + "'.");
					lhs->value.setShared(mVariables[lhs->field]);
				}
				if(rhs->type == VARIABLE) {
					if(!isVariable(rhs->field))
						throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + rhs->field + "'.");
					rhs->value.setShared(mVariables[rhs->field]);
				}
				if(rhs->value.matrix().size() == 1) {
					if(lhs->value.isLocal()) {
						if(op->field == "*" || op->field == ".*")
							lhs->value.local().array() *= rhs->value.matrix()(0, 0);
						else // if(op->field == "/" || op->field == "./")
							lhs->value.local().array() /= rhs->value.matrix()(0, 0);
					} else {
						if(op->field == "*" || op->field == ".*")
							lhs->value.local() = lhs->value.matrix().array() * rhs->value.matrix()(0, 0);
						else // if(op->field == "/" || op->field == "./")
							lhs->value.local() = lhs->value.matrix().array() / rhs->value.matrix()(0, 0);
						lhs->value.mapLocal();
						lhs->type = VALUE;
					}
				} else if(lhs->value.matrix().size() == 1) {
					typename Derived::Scalar temp = lhs->value.matrix()(0, 0);
					if(op->field == "*" || op->field == ".*")
						lhs->value.local() = rhs->value.matrix().array() * temp;
					else // if(op->field == "/" || op->field == "./")
						lhs->value.local() = Derived::Constant(rhs->value.matrix().rows(), rhs->value.matrix().cols(), temp).array() / rhs->value.matrix().array();
					lhs->value.mapLocal();
					lhs->type = VALUE;
				} else if((op->field == ".*" || op->field == "./") && lhs->value.matrix().rows() == rhs->value.matrix().rows() && lhs->value.matrix().cols() == rhs->value.matrix().cols()) {
					if(lhs->value.isLocal()) {
						if(op->field == ".*")
							lhs->value.local().array() *= rhs->value.matrix().array();
						else // if(op->field == "./")
							lhs->value.local().array() /= rhs->value.matrix().array();
					} else {
						if(op->field == ".*")
							lhs->value.local() = lhs->value.matrix().array() * rhs->value.matrix().array();
						else // if(op->field == "./")
							lhs->value.local() = lhs->value.matrix().array() / rhs->value.matrix().array();
						lhs->value.mapLocal();
						lhs->type = VALUE;
					}
				} else if(op->field == "*" && lhs->value.matrix().cols() == rhs->value.matrix().rows()) {
					if(lhs->value.isLocal()) {
						lhs->value.local() *= rhs->value.matrix();
						lhs->value.mapLocal();
					} else {
						lhs->value.local() = lhs->value.matrix() * rhs->value.matrix();
						lhs->value.mapLocal();
						lhs->type = VALUE;
					}
				} else {
					throw std::runtime_error("Invalid operand dimensions for operation '" + lhs->field + op->field + rhs->field + "'.");
				}
				chunks.erase(op, rhs + 1);
				op = lhs + 1;
				rhs = (op != chunks.end()) ? op + 1 : op;
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
				operationPerformed = true;
#	endif
#endif
			} else {
				lhs = op;
				op = rhs;
				rhs++;
			}
		}
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		if(operationPerformed) { std::cout << "*: "; printChunks(chunks); std::cout << std::endl; }
#	endif
#endif
	}

	template <typename Derived>
	void Parser<Derived>::evalAddition(ChunkArray & chunks)
	{
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		bool operationPerformed = false;
#	endif
#endif
		if(chunks.size() < 3) return;
		typename ChunkArray::iterator lhs = chunks.begin(), op = lhs + 1, rhs = op + 1;
		for(; lhs != chunks.end() && op != chunks.end() && rhs != chunks.end();) {
			if((op->type == OPERATOR) && (op->field == "+" || op->field == "-" || op->field == ".+" || op->field == ".-")) {
				if(lhs->type == VARIABLE) {
					if(!isVariable(lhs->field))
						throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + lhs->field + "'.");
					lhs->value.setShared(mVariables[lhs->field]);
				}
				if(rhs->type == VARIABLE) {
					if(!isVariable(rhs->field))
						throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + rhs->field + "'.");
					rhs->value.setShared(mVariables[rhs->field]);
				}
				if(rhs->value.matrix().size() == 1) {
					if(lhs->value.isLocal()) {
						if(op->field == "+" || op->field == ".+")
							lhs->value.local().array() += rhs->value.matrix()(0, 0);
						else // if(op->field == "-" || op->field == ".-")
							lhs->value.local().array() -= rhs->value.matrix()(0, 0);
					} else {
						if(op->field == "+" || op->field == ".+")
							lhs->value.local() = lhs->value.matrix().array() + rhs->value.matrix()(0, 0);
						else // if(op->field == "-" || op->field == ".-")
							lhs->value.local() = lhs->value.matrix().array() - rhs->value.matrix()(0, 0);
						lhs->value.mapLocal();
						lhs->type = VALUE;
					}
				} else if(lhs->value.matrix().size() == 1) {
					typename Derived::Scalar temp = lhs->value.matrix()(0, 0);
					if(op->field == "+" || op->field == ".+")
						lhs->value.local() = rhs->value.matrix().array() + temp;
					else // if(op->field == "-" || op->field == ".-")
						lhs->value.local() = Derived::Constant(rhs->value.matrix().rows(), rhs->value.matrix().cols(), temp).array() - rhs->value.matrix().array();
					lhs->value.mapLocal();
					lhs->type = VALUE;
				} else if(lhs->value.matrix().rows() == rhs->value.matrix().rows() && lhs->value.matrix().cols() == rhs->value.matrix().cols()) {
					if(lhs->value.isLocal()) {
						if(op->field == "+" || op->field == ".+")
							lhs->value.local().array() += rhs->value.matrix().array();
						else // if(op->field == "-" || op->field == ".-")
							lhs->value.local().array() -= rhs->value.matrix().array();
					} else {
						if(op->field == "+" || op->field == ".+")
							lhs->value.local() = lhs->value.matrix().array() + rhs->value.matrix().array();
						else // if(op->field == "-" || op->field == ".-")
							lhs->value.local() = lhs->value.matrix().array() - rhs->value.matrix().array();
						lhs->value.mapLocal();
						lhs->type = VALUE;
					}
				} else {
					throw std::runtime_error("Invalid operand dimensions for operation '" + lhs->field + op->field + rhs->field + "'.");
				}
				chunks.erase(op, rhs + 1);
				op = lhs + 1;
				rhs = (op != chunks.end()) ? op + 1 : op;
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
				operationPerformed = true;
#	endif
#endif
			} else {
				lhs = op;
				op = rhs;
				rhs++;
			}
		}
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		if(operationPerformed) { std::cout << "+: "; printChunks(chunks); std::cout << std::endl; }
#	endif
#endif
	}

	template <typename Derived>
	void Parser<Derived>::evalAssignment(ChunkArray & chunks)
	{
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		bool operationPerformed = false;
#	endif
#endif
		if(chunks.size() < 3) return;
		typename ChunkArray::iterator rhs = chunks.end() - 1, op = rhs - 1, lhs = op - 1;
		for(; op != chunks.begin() && rhs != chunks.begin();) {
			if(op->type == OPERATOR && op->field == "=" && (lhs->type == VALUE || lhs->type == VARIABLE) && (rhs->type == VALUE || rhs->type == VARIABLE)) {
				if(rhs->type == VARIABLE) {
					if(!isVariable(rhs->field))
						throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + rhs->field + "'.");
					rhs->value.setShared(mVariables[rhs->field]);
				}
				if(lhs->type == VALUE) {
					lhs->value.local() = rhs->value.matrix();
					lhs->value.mapLocal();
				} else { //if(lhs->type == VARIABLE) {
					if(isVariable(lhs->field)) {
						lhs->value.setShared(mVariables[lhs->field]);
						if(lhs->row0 == -1) {
							if(lhs->value.matrix().rows() == rhs->value.matrix().rows() && lhs->value.matrix().cols() == rhs->value.matrix().cols()) {
								lhs->value.matrix() = rhs->value.matrix();
							} else {
								mVariables[lhs->field].local() = rhs->value.matrix();
								mVariables[lhs->field].mapLocal();
							}
						} else { //if(lhs->row0 != -1) {
							lhs->value.matrix().block(lhs->row0, lhs->col0, lhs->rows, lhs->cols) = rhs->value.matrix();
						}
					} else {
						mVariables[lhs->field].local() = rhs->value.matrix();
						mVariables[lhs->field].mapLocal();
					}
				}
				rhs = chunks.erase(op, rhs + 1);
				op = (rhs != chunks.begin()) ? rhs - 1 : rhs;
				if (op != chunks.begin()) lhs = op - 1;
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
				operationPerformed = true;
#	endif
#endif
			} else {
				rhs = op;
				op = lhs;
				lhs--;
			}
		}
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
		if(operationPerformed) { std::cout << "=: "; printChunks(chunks); std::cout << std::endl; }
#	endif
#endif
	}
	
#ifdef DEBUG
#	ifdef EIGENLAB_DEBUG
	template <typename Derived>
	void Parser<Derived>::printChunks(ChunkArray & chunks, size_t maxRows, size_t maxCols, int precision)
	{
		std::cout << "__";
		for(typename ChunkArray::iterator it = chunks.begin(); it != chunks.end(); it++) {
			switch(it->type) {
				case VALUE:
					std::cout << textRepresentation(it->value, maxRows, maxCols, precision);
					if(it->row0 != -1)
						std::cout << "(" << it->row0 << ":" << it->row0 + it->rows - 1 << "," << it->col0 << ":" << it->col0 + it->cols - 1 << ")";
					break;
				case VARIABLE:
					std::cout << it->field;// << "=" << textRepresentation(mVariables[it->field], maxRows, maxCols, precision);
					if(it->row0 != -1)
						std::cout << "(" << it->row0 << ":" << it->row0 + it->rows - 1 << "," << it->col0 << ":" << it->col0 + it->cols - 1 << ")";
					break;
				case OPERATOR:
					std::cout << it->field;
					break;
				case FUNCTION:
					std::cout << "f()=" << it->field;
					break;
			}
			std::cout << "__";
		}
	}
	
	template <typename Derived>
	void Parser<Derived>::printVars(size_t maxRows, size_t maxCols, int precision)
	{
		for(typename ValueMap::iterator it = mVariables.begin(); it != mVariables.end(); it++)
			std::cout << it->first << " (" << it->second.matrix().rows() << "x" << it->second.matrix().cols() << ") = " << textRepresentation(it->second, maxRows, maxCols, precision) << std::endl;
	}
	
	template <typename Derived>
	std::string Parser<Derived>::textRepresentation(Value<Derived> & val, size_t maxRows, size_t maxCols, int precision)
	{
		if(val.matrix().size() == 1)
			return numberToString<typename Derived::Scalar>(val.matrix()(0, 0), precision);
		else {
			std::string str = "[";
			for(size_t row = 0; row < val.matrix().rows() && row < maxRows; row++) {
				str += (row > 0 ? ";[" : "[");
				for(size_t col = 0; col < val.matrix().cols() && col < maxCols; col++) {
					if(col > 0) str += ",";
					str += numberToString<typename Derived::Scalar>(val.matrix()(row, col), precision);
				}
				str += (val.matrix().cols() > maxCols ? "...]" : "]");
			}
			str += (val.matrix().rows() > maxRows ? "...]" : "]");
			return str;
		}
	}
#	endif // #ifdef EIGENLAB_DEBUG
#endif // #ifdef DEBUG
	
	template <typename Derived>
	std::string Parser<Derived>::trim(const std::string & str)
	{
		if(str.size() == 0) return str;
		std::string::const_iterator first = str.begin();
		std::string::const_iterator last = str.end() - 1;
		while(first < last && isspace(*first)) first++;
		while(first < last && isspace(*last)) last--;
		return std::string(first, last + 1);
	}
	
	template <typename Derived>
	std::vector<std::string> Parser<Derived>::split(const std::string & str, const char delimeter)
	{
		std::vector<std::string> args;
		std::string::const_iterator i0 = str.begin();
		for(std::string::const_iterator it = str.begin(); it != str.end(); it++) {
			if((*it) == delimeter) {
				args.push_back(trim(std::string(i0, it)));
				i0 = it + 1;
			}
		}
		args.push_back(std::string(i0, str.end()));
		return args;
	}
	
	template <typename Derived>
	template <typename T>
	bool Parser<Derived>::isNumber(const std::string & str, T * num)
	{
		std::istringstream iss(str);
		if(num)
			iss >> (* num);
		else {
			T number;
			iss >> number;
		}
		return (!iss.fail() && !iss.bad() && iss.eof());
	}
	
	template <typename Derived>
	template<typename T>
	T Parser<Derived>::stringToNumber(const std::string & str)
	{
		std::istringstream iss(str);
		T number;
		iss >> number;
		if(iss.fail() || iss.bad() || !iss.eof())
			throw std::runtime_error("Failed to convert " + str + " to a number.");
		return number;
	}
	
	template <typename Derived>
	template<typename T>
	std::string Parser<Derived>::numberToString(T num, int precision)
	{
		std::ostringstream oss;
		if(precision)
			oss << std::setprecision(precision) << num;
		else
			oss << num;
		return oss.str();
	}
	
#ifdef DEBUG
	template <typename Derived>
	void
	Parser<Derived>::test_w_lt(size_t & numFails,
				   typename Derived::Scalar & /* s */,
				   Derived & a34,
				   Derived & /* b34 */,
				   Derived & /* c43 */,
				   Derived & /* v */, std::true_type)
	{
		//
		// tests that only work if Derived::Scalar has operator<
		//
		Value<Derived> resultValue;
		Derived resultMatrix;
		Derived temp;
		std::cout << "Test min(a): ";
		resultValue = eval("min(a)");
		resultMatrix.setConstant(1, 1, a34.minCoeff());
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }

		std::cout << "Test min(a, 0): ";
		resultValue = eval("min(a, 0)");
		resultMatrix = a34.colwise().minCoeff();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }

		std::cout << "Test min(a, 1): ";
		resultValue = eval("min(a, 1)");
		resultMatrix = a34.rowwise().minCoeff();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }

		std::cout << "Test max(a): ";
		resultValue = eval("max(a)");
		resultMatrix.setConstant(1, 1, a34.maxCoeff());
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }

		std::cout << "Test max(a, 0): ";
		resultValue = eval("max(a, 0)");
		resultMatrix = a34.colwise().maxCoeff();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }

		std::cout << "Test max(a, 1): ";
		resultValue = eval("max(a, 1)");
		resultMatrix = a34.rowwise().maxCoeff();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }

		std::cout << "Test absmax(a): ";
		resultValue = eval("absmax(a)");
		resultMatrix.setConstant(1, 1, std::abs(a34.maxCoeff()) >= std::abs(a34.minCoeff()) ? a34.maxCoeff() : a34.minCoeff());
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }

		std::cout << "Test absmax(a, 0): ";
		resultValue = eval("absmax(a, 0)");
		resultMatrix = a34.colwise().maxCoeff();
		temp = a34.colwise().minCoeff();
		for(Index i = 0; i < resultMatrix.size(); ++i) {
			if(std::abs(resultMatrix(i)) < std::abs(temp(i)))
				resultMatrix(i) = temp(i);
		}
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }

		std::cout << "Test absmax(a, 1): ";
		resultValue = eval("absmax(a, 1)");
		resultMatrix = a34.rowwise().maxCoeff();
		temp = a34.rowwise().minCoeff();
		for(Index i = 0; i < resultMatrix.size(); ++i) {
			if(std::abs(resultMatrix(i)) < std::abs(temp(i)))
				resultMatrix(i) = temp(i);
		}
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
	}

	template <typename Derived>
	void
	Parser<Derived>::test_w_lt(size_t & /* numFails */,
				   typename Derived::Scalar & /* s */,
				   Derived & /* a34 */,
				   Derived & /* b34 */,
				   Derived & /* c43 */,
				   Derived & /* v */, std::false_type)
	{
		// do nothing
	}

	template <typename Derived>
	size_t Parser<Derived>::test()
	{
		std::cout << std::endl;
		std::cout << "BEGIN unit test for EigenLab..." << std::endl;
		std::cout << "Make sure this function completes successfuly and prints the message "
			"'Successfully completed unit test for EigenLab with no failures.'" << std::endl;
		std::cout << std::endl;
		
		size_t numFails = 0;
		Value<Derived> resultValue;
		Derived resultMatrix;
		Derived temp;
		typename Derived::Scalar s = 2;
		
		Derived a34 = Derived::Random(3, 4);
		Derived b34 = Derived::Random(3, 4);
		Derived c43 = Derived::Random(4, 3);
		Derived v = Derived::Random(1, 10);
		//std::cout << "a34=" << std::endl << a34 << std::endl << std::endl;
		//std::cout << "b34=" << std::endl << b34 << std::endl << std::endl;
		//std::cout << "c43=" << std::endl << c43 << std::endl << std::endl;
		//std::cout << "v=" << std::endl << v << std::endl << std::endl;
		
		var("a").setShared(a34);
		var("b").setShared(b34);
		var("c").setShared(c43);
		var("v").setShared(v);
		var("s").setShared(& s);
		
		////////////////////////////////////////
		std::cout << "Testing basic operations..." << std::endl << std::endl;
		////////////////////////////////////////
		
		std::cout << "Test matrix addition a + b: ";
		resultValue = eval("a + b");
		resultMatrix = a34 + b34;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix/scalar addition a + s: ";
		resultValue = eval("a + s");
		resultMatrix = a34.array() + s;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test scalar/matrix addition s + b: ";
		resultValue = eval("s + b");
		resultMatrix = b34.array() + s;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix addition a .+ b: ";
		resultValue = eval("a .+ b");
		resultMatrix = a34 + b34;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix/scalar addition a .+ s: ";
		resultValue = eval("a .+ s");
		resultMatrix = a34.array() + s;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test scalar/matrix addition s .+ b: ";
		resultValue = eval("s .+ b");
		resultMatrix = b34.array() + s;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix subtraction a - b: ";
		resultValue = eval("a - b");
		resultMatrix = a34 - b34;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix/scalar subtraction a - s: ";
		resultValue = eval("a - s");
		resultMatrix = a34.array() - s;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test scalar/matrix subtraction s - b: ";
		resultValue = eval("s - b");
		resultMatrix = (-b34.array()) + s;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix subtraction a .- b: ";
		resultValue = eval("a .- b");
		resultMatrix = a34 - b34;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix/scalar subtraction a .- s: ";
		resultValue = eval("a .- s");
		resultMatrix = a34.array() - s;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test scalar/matrix subtraction s .- b: ";
		resultValue = eval("s .- b");
		resultMatrix = (-b34.array()) + s;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix negation -a: ";
		resultValue = eval("-a");
		resultMatrix = -a34;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test scalar negation -s: ";
		resultValue = eval("-s");
		resultMatrix.setConstant(1, 1, -s);
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix coefficient-wise multiplication a .* b: ";
		resultValue = eval("a .* b");
		resultMatrix = a34.array() * b34.array();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix/scalar coefficient-wise multiplication a * s: ";
		resultValue = eval("a * s");
		resultMatrix = a34.array() * s;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test scalar/matrix coefficient-wise multiplication s * b: ";
		resultValue = eval("s * b");
		resultMatrix = b34.array() * s;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix/scalar coefficient-wise multiplication a .* s: ";
		resultValue = eval("a .* s");
		resultMatrix = a34.array() * s;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test scalar/matrix coefficient-wise multiplication s .* b: ";
		resultValue = eval("s .* b");
		resultMatrix = b34.array() * s;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix coefficient-wise division a ./ b: ";
		resultValue = eval("a ./ b");
		resultMatrix = a34.array() / b34.array();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix/scalar coefficient-wise division a / s: ";
		resultValue = eval("a / s");
		resultMatrix = a34.array() / s;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test scalar/matrix coefficient-wise division s / b: ";
		resultValue = eval("s / b");
		resultMatrix = Derived::Constant(b34.rows(), b34.cols(), s).array() / b34.array();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix/scalar coefficient-wise division a ./ s: ";
		resultValue = eval("a ./ s");
		resultMatrix = a34.array() / s;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test scalar/matrix coefficient-wise division s ./ b: ";
		resultValue = eval("s ./ b");
		resultMatrix = Derived::Constant(b34.rows(), b34.cols(), s).array() / b34.array();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix coefficient-wise power a .^ b: ";
		resultValue = eval("abs(a) .^ b");
		resultMatrix = a34;
		for(Index i = 0; i < a34.size(); ++i)
			resultMatrix(i) = pow(std::abs(a34(i)), b34(i));
		//		std::cout << std::endl;
		//		std::cout << "a=" << std::endl << a34 << std::endl << std::endl;
		//		std::cout << "b=" << std::endl << b34 << std::endl << std::endl;
		//		std::cout << "val=" << std::endl << resultValue.matrix() << std::endl << std::endl;
		//		std::cout << "mat=" << std::endl << resultMatrix << std::endl << std::endl;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix/scalar coefficient-wise power a ^ s: ";
		resultValue = eval("abs(a) ^ s");
		resultMatrix = a34;
		for(Index i = 0; i < a34.size(); ++i)
			resultMatrix(i) = pow(std::abs(a34(i)), s);
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test scalar/matrix coefficient-wise power s ^ b: ";
		resultValue = eval("s ^ b");
		resultMatrix = b34;
		for(Index i = 0; i < b34.size(); ++i)
			resultMatrix(i) = pow(s, b34(i));
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix/scalar coefficient-wise power a .^ s: ";
		resultValue = eval("abs(a) .^ s");
		resultMatrix = a34;
		for(Index i = 0; i < a34.size(); ++i)
			resultMatrix(i) = pow(std::abs(a34(i)), s);
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test scalar/matrix coefficient-wise power s .^ b: ";
		resultValue = eval("s .^ b");
		resultMatrix = b34;
		for(Index i = 0; i < b34.size(); ++i)
			resultMatrix(i) = pow(s, b34(i));
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix multiplication a * b: ";
		resultValue = eval("a * c");
		resultMatrix = a34 * c43;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test grouped subexpression (a + b) * c: ";
		resultValue = eval("(a + b) * c");
		resultMatrix = (a34 + b34) * c43;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test nested groups ((a + (b * 3 + 1)) * c).^2: ";
		resultValue = eval("((a + (b * 3 + 1)) * c).^2");
		resultMatrix = ((a34.array() + (b34.array() * 3 + 1)).matrix() * c43).array().pow(2);
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		////////////////////////////////////////
		std::cout << std::endl << "Testing coefficient and submatrix block access..." << std::endl << std::endl;
		////////////////////////////////////////
		
		std::cout << "Test matrix coefficient access a(i,j): ";
		resultValue = eval("a(1,2)");
		resultMatrix.setConstant(1, 1, a34(1,2));
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test submatrix block access a(i:p,j:q): ";
		resultValue = eval("a(1:2,2:3)");
		resultMatrix = a34.block(1, 2, 2, 2);
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test submatrix block access using 'end' and ':' identifiers a(i:end,:): ";
		resultValue = eval("a(1:end,:)");
		resultMatrix = a34.block(1, 0, a34.rows() - 1, a34.cols());
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test submatrix block access using subexpressions: ";
		resultValue = eval("a(2-1:2-1,0+1:3-1)");
		resultMatrix = a34.block(1, 1, 1, 2);
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test submatrix block access using subexpressions with 'end' keyword: ";
		resultValue = eval("a(2-1:end-1,0+1:end-1)");
		resultMatrix = a34.block(1, 1, a34.rows() - 2, a34.cols() - 2);
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test vector coefficient access v(i): ";
		resultValue = eval("v(5)");
		resultMatrix.setConstant(1, 1, v(5));
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test subvector segment access v(i:j): ";
		resultValue = eval("v(3:6)");
		resultMatrix = v.block(0, 3, 1, 4);
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test subvector segment access using 'end' identifier v(i:end): ";
		resultValue = eval("v(5:end)");
		resultMatrix = v.block(0, 5, 1, v.cols() - 5);
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test subvector segment access using ':' identifier v(:): ";
		resultValue = eval("v(:)");
		resultMatrix = v;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test subvector segment access using subexpressions: ";
		resultValue = eval("v(3-1:5+2)");
		resultMatrix = v.block(0, 2, 1, 6);
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test subvector segment access using subexpressions with 'end' keyword: ";
		resultValue = eval("v((end-8)*2:end-3)");
		resultMatrix = v.block(0, 2, 1, 5);
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		////////////////////////////////////////
		std::cout << std::endl << "Testing vector/matrix expressions..." << std::endl << std::endl;
		////////////////////////////////////////
		
		std::cout << "Test numeric range expression [i:j]: ";
		resultValue = eval("[2:5]");
		resultMatrix.resize(1, 4);
		resultMatrix << 2, 3, 4, 5;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test numeric range expression [i:s:j]: ";
		resultValue = eval("[2:2:10]");
		resultMatrix.resize(1, 5);
		resultMatrix << 2, 4, 6, 8, 10;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test numeric range with subexpressions: ";
		resultValue = eval("[6-2:5*2-3]");
		std::cout << "val=" << std::endl << resultValue.matrix() << std::endl << std::endl;
		resultMatrix.resize(1, 4);
		resultMatrix << 4, 5, 6, 7;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix expression [[1, 2]; [3, 4]]: ";
		resultValue = eval("[[1, 2]; [3, 4]]");
		resultMatrix.resize(2, 2);
		resultMatrix << 1, 2, 3, 4;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test matrix expression [  1, 2, 3;	4:6 ]: ";
		resultValue = eval("[1, 2, 3;	4:6]");
		resultMatrix.resize(2, 3);
		resultMatrix << 1, 2, 3, 4, 5, 6;
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		////////////////////////////////////////
		std::cout << std::endl << "Testing coefficient-wise functions..." << std::endl << std::endl;
		////////////////////////////////////////
		
		std::cout << "Test coefficient-wise abs(a): ";
		resultValue = eval("abs(a)");
		resultMatrix.resize(3, 4);
		resultMatrix.real() = a34.array().abs();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test coefficient-wise sqrt(a): ";
		resultValue = eval("sqrt(abs(a))");
		resultMatrix.real() = a34.array().abs().sqrt();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test coefficient-wise exp(a): ";
		resultValue = eval("exp(abs(a) + 0.001)");
		resultMatrix.real() = (a34.array().abs() + 0.001).exp();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test coefficient-wise log(a): ";
		resultValue = eval("log(abs(a) + 0.001)");
		resultMatrix.real() = (a34.array().abs() + 0.001).log();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test coefficient-wise log10(a): ";
		resultValue = eval("log10(abs(a) + 0.001)");
		resultMatrix.real() = (a34.array().abs() + 0.001).log() * (1.0 / log(10));
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test coefficient-wise sin(a): ";
		resultValue = eval("sin(a)");
		resultMatrix = a34.array().sin();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test coefficient-wise cos(a): ";
		resultValue = eval("cos(a)");
		resultMatrix = a34.array().cos();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test coefficient-wise tan(a): ";
		resultValue = eval("tan(a)");
		resultMatrix = a34.array().tan();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test coefficient-wise asin(a): ";
		resultValue = eval("asin(a)");
		resultMatrix = a34.array().asin();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test coefficient-wise acos(a): ";
		resultValue = eval("acos(a)");
		resultMatrix = a34.array().acos();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		////////////////////////////////////////
		std::cout << std::endl << "Testing matrix reduction functions..." << std::endl << std::endl;
		////////////////////////////////////////
		
		std::cout << "Test trace(a): ";
		resultValue = eval("trace(a)");
		resultMatrix.setConstant(1, 1, a34.trace());
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test norm(a): ";
		resultValue = eval("norm(a)");
		resultMatrix.setConstant(1, 1, a34.norm());
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test size(a, 0): ";
		resultValue = eval("size(a, 0)");
		resultMatrix.setConstant(1, 1, a34.rows());
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test size(a, 1): ";
		resultValue = eval("size(a, 1)");
		resultMatrix.setConstant(1, 1, a34.cols());
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		test_w_lt(numFails, s, a34, b34, c43, v, has_operator_lt<typename Derived::Scalar>());
		
		std::cout << "Test mean(a): ";
		resultValue = eval("mean(a)");
		resultMatrix.setConstant(1, 1, a34.mean());
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test mean(a, 0): ";
		resultValue = eval("mean(a, 0)");
		resultMatrix = a34.colwise().mean();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test mean(a, 1): ";
		resultValue = eval("mean(a, 1)");
		resultMatrix = a34.rowwise().mean();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test sum(a): ";
		resultValue = eval("sum(a)");
		resultMatrix.setConstant(1, 1, a34.sum());
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test sum(a, 0): ";
		resultValue = eval("sum(a, 0)");
		resultMatrix = a34.colwise().sum();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test sum(a, 1): ";
		resultValue = eval("sum(a, 1)");
		resultMatrix = a34.rowwise().sum();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test prod(a): ";
		resultValue = eval("prod(a)");
		resultMatrix.setConstant(1, 1, a34.prod());
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test prod(a, 0): ";
		resultValue = eval("prod(a, 0)");
		resultMatrix = a34.colwise().prod();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test prod(a, 1): ";
		resultValue = eval("prod(a, 1)");
		resultMatrix = a34.rowwise().prod();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		////////////////////////////////////////
		std::cout << std::endl << "Testing matrix functions..." << std::endl << std::endl;
		////////////////////////////////////////
		
		std::cout << "Test transpose(a): ";
		resultValue = eval("transpose(a)");
		resultMatrix = a34.transpose();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test conjugate(a): ";
		resultValue = eval("conjugate(a)");
		resultMatrix = a34.conjugate();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test adjoint(a): ";
		resultValue = eval("adjoint(a)");
		resultMatrix = a34.adjoint();
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
        else { std::cout << "FAIL" << std::endl; ++numFails; }
        
        ////////////////////////////////////////
        std::cout << std::endl << "Testing matrix initializers..." << std::endl << std::endl;
        ////////////////////////////////////////
        
        std::cout << "Test zeros(5): ";
        resultValue = eval("zeros(5)");
        resultMatrix = Derived::Zero(5, 5);
        if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
        else { std::cout << "FAIL" << std::endl; ++numFails; }
        
        std::cout << "Test ones(5): ";
        resultValue = eval("ones(5)");
        resultMatrix = Derived::Ones(5, 5);
        if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
        else { std::cout << "FAIL" << std::endl; ++numFails; }
        
        std::cout << "Test eye(5): ";
        resultValue = eval("eye(5)");
        resultMatrix = Derived::Identity(5, 5);
        if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
        else { std::cout << "FAIL" << std::endl; ++numFails; }
        
        try {
            std::cout << "Test zeros(5.2): ";
            resultValue = eval("zeros(5.2)"); // <-- Should NOT succeed!!!
            std::cout << "FAIL" << std::endl; ++numFails;
        } catch(std::runtime_error &err) {
            std::cout << err.what() << std::endl;
            std::cout << "Exception caught, so we're OK" << std::endl;
        }
        
        try {
            std::cout << "Test eye(-3): ";
            resultValue = eval("eye(-3)"); // <-- Should NOT succeed!!!
            std::cout << "FAIL" << std::endl; ++numFails;
        } catch(std::runtime_error &err) {
            std::cout << err.what() << std::endl;
            std::cout << "Exception caught, so we're OK" << std::endl;
        }
        
        std::cout << "Test zeros(4,7): ";
        resultValue = eval("zeros(4,7)");
        resultMatrix = Derived::Zero(4, 7);
        if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
        else { std::cout << "FAIL" << std::endl; ++numFails; }
        
        std::cout << "Test ones(4,7): ";
        resultValue = eval("ones(4,7)");
        resultMatrix = Derived::Ones(4, 7);
        if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
        else { std::cout << "FAIL" << std::endl; ++numFails; }
        
        std::cout << "Test eye(4,7): ";
        resultValue = eval("eye(4,7)");
        resultMatrix = Derived::Identity(4, 7);
        if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
        else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		////////////////////////////////////////
		std::cout << std::endl << "Testing variable assignment..." << std::endl << std::endl;
		////////////////////////////////////////
		
		std::cout << "Test assigning to a variable with the same dimensions a = b: ";
		eval("a = b");
		if(a34.isApprox(b34)) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test assigning to a variable with different dimensions a = c: ";
		eval("a = c");
		if(var("a").matrix().isApprox(c43) && a34.isApprox(b34)) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		var("a").setShared(a34);
		
		std::cout << "Test creating a new variable x = [1,2;3,4]: ";
		resultValue = eval("x = [1,2;3,4]");
		if(var("x").matrix().isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test assigning to a variable coefficient a(i,j) = s: ";
		eval("a(1, 2) = s");
		if(a34(1, 2) == s) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
		std::cout << "Test assigning to a variable submatrix block a(0:1,1:2) = x: ";
		resultValue = eval("a(0:1,1:2) = x");
		if(a34.block(0,1,2,2).isApprox(var("x").matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }
		
        try {
            std::cout << "Test bad function call: ";
            resultValue = eval("foobar(-3)"); // <-- Should NOT succeed!!!
            std::cout << "FAIL" << std::endl; ++numFails;
        } catch(std::runtime_error &err) {
            std::cout << err.what() << std::endl;
            std::cout << "Exception caught, so we're OK" << std::endl;
        }
		////////////////////////////////////////
		std::cout << std::endl << "Testing fp parsing..." << std::endl << std::endl;
		////////////////////////////////////////

		std::cout << "Test assigning 1.2e-3: ";
		resultValue = eval("s = 1.2e-3");
		resultMatrix.setConstant(1, 1, typename Derived::Scalar(1.2e-3));
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }

		std::cout << "Test assigning 1.e3: ";
		resultValue = eval("s = 1.e3");
		resultMatrix.setConstant(1, 1, typename Derived::Scalar(1000));
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }

		std::cout << "Test assigning 12.34e05: ";
		resultValue = eval("s = 12.34e05");
		resultMatrix.setConstant(1, 1, typename Derived::Scalar(123.4e4));
		if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
		else { std::cout << "FAIL" << std::endl; ++numFails; }

		std::cout << std::endl;
		if(numFails == 0)
			std::cout << "Successfully completed unit test for EigenLab with no failures." << std::endl;
		else
			std::cout << "Completed unit test for EigenLab with " << numFails << " failures (see above)." << std::endl;
		std::cout << std::endl;
		return numFails;
	}
#endif // #ifdef DEBUG

} // namespace EigenLab

#endif // #ifndef EigenLab_H