.. _program_listing_file__tmp_ws_src_grid_map_grid_map_filters_include_EigenLab_EigenLab.hpp:

Program Listing for File EigenLab.hpp
=====================================

|exhale_lsh| :ref:`Return to documentation for file <file__tmp_ws_src_grid_map_grid_map_filters_include_EigenLab_EigenLab.hpp>` (``/tmp/ws/src/grid_map/grid_map_filters/include/EigenLab/EigenLab.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   // --*-  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__EIGENLAB_HPP_
   #define EIGENLAB__EIGENLAB_HPP_
   
   #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)
     {
     }
     explicit Value(const typename Derived::Scalar s)
     : mLocal(Derived::Constant(1, 1, s)), mShared(
         mLocal.data(), mLocal.rows(), mLocal.cols()), mIsLocal(true)
     {
     }
     explicit 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.
     explicit 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>
     // deepcode ignore CopyPasteError: We will not change third party code yet.
     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("square");
     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("minOfFinites");
       mFunctions.push_back("max");
       mFunctions.push_back("maxOfFinites");
       mFunctions.push_back("absmax");
       mFunctions.push_back("cwiseMin");
       mFunctions.push_back("cwiseMax");
     }
     mFunctions.push_back("mean");
     mFunctions.push_back("meanOfFinites");
     mFunctions.push_back("sum");
     mFunctions.push_back("sumOfFinites");
     mFunctions.push_back("prod");
     mFunctions.push_back("numberOfFinites");
   
     // 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 = static_cast<int>(chunks.back().value.matrix().rows());
         int cols = static_cast<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;
             // (last == -1 ? int(chunks.back().value.matrix().rows()) : last + 1) - first;
             chunks.back().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;
             // (last == -1 ? int(chunks.back().value.matrix().cols()) : last + 1) - first;
             chunks.back().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;
           // (last == -1 ? int(chunks.back().value.matrix().rows()) : last + 1) - first;
           chunks.back().rows = last + 1 - first;
           evalIndexRange(args[1], &first, &last, cols);
           chunks.back().col0 = first;
           // (last == -1 ? int(chunks.back().value.matrix().cols()) : last + 1) - first;
           chunks.back().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;}
           static const 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, static_cast<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[static_cast<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[static_cast<int>(row0) - 1].size()) {
         throw std::runtime_error(
                 "Invalid matrix definition '[" + str + "]'. Successive row entries '" +
                 rows[static_cast<int>(i) - 1] + "' and '" + rows[i] +
                 "' do not have the same number of columns.");
       }
       row0 += nrows;
     }
     if (temp.empty()) {return;}
     nrows = temp.size();
     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 == "minOfFinites") {
       result.setLocal(arg.matrix().minCoeffOfFinites());
       return true;
     } else if (name == "max") {
       result.setLocal(arg.matrix().maxCoeff());
       return true;
     } else if (name == "maxOfFinites") {
       result.setLocal(arg.matrix().maxCoeffOfFinites());
       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;
       }
     } else if (name == "cwiseMin") {
       if (arg1.matrix().size() == 1) {
         typename Derived::Scalar arg1scalar = arg1.matrix()(0, 0);
         Derived arg1matrix = Derived::Constant(
           arg0.matrix().rows(),
           arg0.matrix().cols(), arg1scalar);
         result.local() = arg0.matrix().cwiseMin(arg1matrix);
         result.mapLocal();
         return true;
       } else if (  // NOLINT
         arg0.matrix().cols() == arg1.matrix().cols() &&
         arg0.matrix().rows() == arg1.matrix().rows())
       {
         result.local() = arg0.matrix().cwiseMin(arg1.matrix());
         result.mapLocal();
         return true;
       } else {
         throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
       }
     } else if (name == "cwiseMax") {
       if (arg1.matrix().size() == 1) {
         typename Derived::Scalar arg1scalar = arg1.matrix()(0, 0);
         Derived arg1matrix = Derived::Constant(
           arg0.matrix().rows(),
           arg0.matrix().cols(), arg1scalar);
         result.local() = arg0.matrix().cwiseMax(arg1matrix);
         result.mapLocal();
         return true;
       } else if (  // NOLINT
         arg0.matrix().cols() == arg1.matrix().cols() &&
         arg0.matrix().rows() == arg1.matrix().rows())
       {
         result.local() = arg0.matrix().cwiseMax(arg1.matrix());
         result.mapLocal();
         return true;
       } else {
         throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
       }
     }
     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 == "square") {
         result.local() = arg.matrix().array().square();
         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 (  // NOLINT
         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 == "meanOfFinites") {
         result.setLocal(arg.matrix().meanOfFinites());
         return;
       } else if (name == "sum") {
         result.setLocal(arg.matrix().sum());
         return;
       } else if (name == "sumOfFinites") {
         result.setLocal(arg.matrix().sumOfFinites());
         return;
       } else if (name == "prod") {
         result.setLocal(arg.matrix().prod());
         return;
       } else if (name == "numberOfFinites") {
         result.setLocal(arg.matrix().numberOfFinites());
         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 (  // NOLINT
         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 (  // NOLINT
           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 (  // NOLINT
           (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 (  // NOLINT
           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.empty()) {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;
     }
   
     std::cout << "Test cwiseMin(a, b): ";
     resultValue = eval("cwiseMin(a, b)");
     resultMatrix = a34.cwiseMin(b34);
     if (resultMatrix.isApprox(resultValue.matrix())) {std::cout << "OK" << std::endl;} else {
       std::cout << "FAIL" << std::endl; ++numFails;
     }
   
     std::cout << "Test cwiseMax(a, b): ";
     resultValue = eval("cwiseMax(a, b)");
     resultMatrix = a34.cwiseMax(b34);
     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  // DEBUG
   
   }  // namespace EigenLab
   #endif  // EIGENLAB__EIGENLAB_HPP_