// --*- Mode: C++; c-basic-offset:4; indent-tabs-mode:t; tab-width:4 -*--
// EigenLab
// Version: 1.0.0
// Author: Dr. Marcel Paz Goldschen-Ohm
// Email: marcel.goldschen@gmail.com
// Copyright (c) 2015 by Dr. Marcel Paz Goldschen-Ohm.
// Licence: MIT
//----------------------------------------
#ifndef EigenLab_H
#define EigenLab_H
#include
#include
#include
#include
#include
#include
#include
#include
// 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
#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
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 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 & matrix() { return mShared; }
inline const Eigen::Map & 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(mLocal.data(), mLocal.rows(), mLocal.cols()); mIsLocal = true; }
// Copy shared data to local data (if needed).
inline void copySharedToLocal() { if(!isLocal()) { mLocal = mShared; mapLocal(); } }
// Set local data.
Value() : mLocal(1, 1), mShared(mLocal.data(), mLocal.rows(), mLocal.cols()), mIsLocal(true) {}
Value(const typename Derived::Scalar s) : mLocal(Derived::Constant(1, 1, s)), mShared(mLocal.data(), mLocal.rows(), mLocal.cols()), mIsLocal(true) {}
Value(const Derived & mat) : mLocal(mat), mShared(mLocal.data(), mLocal.rows(), mLocal.cols()), mIsLocal(true) {}
inline void setLocal(const typename Derived::Scalar s) { mLocal = Derived::Constant(1, 1, s); mapLocal(); }
inline void setLocal(const Eigen::MatrixBase & mat) { mLocal = mat; mapLocal(); }
inline void setLocal(const Value & val) { mLocal = val.matrix(); mapLocal(); }
inline void setLocal(const typename Derived::Scalar * data, size_t rows = 1, size_t cols = 1) { setShared(data, rows, cols); copySharedToLocal(); }
inline Value & operator = (const typename Derived::Scalar s) { setLocal(s); return (* this); }
inline Value & operator = (const Derived & mat) { setLocal(mat); return (* this); }
// Set shared data.
Value(const typename Derived::Scalar * data, size_t rows = 1, size_t cols = 1) : mShared(const_cast(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(const_cast(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 ValueXd;
typedef Value ValueXf;
typedef Value ValueXi;
// check if a class has a comparison operator (ie. std::complex does not)
template
struct has_operator_lt_impl
{
template
static auto test(U*) -> decltype(std::declval() < std::declval());
template
static auto test(...) -> std::false_type;
using type = typename std::is_same::type;
};
template
struct has_operator_lt : has_operator_lt_impl::type {};
//----------------------------------------
// Equation parser.
//
// Template typename Derived can be any dynamically sized matrix type supported by Eigen.
//----------------------------------------
template
class Parser
{
public:
// A map to hold named values.
typedef std::map 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 mFunctions;
// Expressions are parsed by first splitting them into chunks.
struct Chunk {
std::string field;
int type;
Value value;
int row0, col0, rows, cols;
Chunk(const std::string & str = "", int t = -1, const Value & val = Value()) : field(str), type(t), value(val), row0(-1), col0(-1), rows(-1), cols(-1) {}
};
enum ChunkType { VALUE = 0, VARIABLE, OPERATOR, FUNCTION };
typedef std::vector ChunkArray;
typedef typename Derived::Index Index;
bool mCacheChunkedExpressions;
std::map 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 & 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 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 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 & mat);
void evalFunction(const std::string & name, std::vector & args, Value & result);
bool evalFunction_1_lt(const std::string & name, Value & arg, Value & result, std::false_type);
bool evalFunction_1_lt(const std::string & name, Value & arg, Value & result, std::true_type);
bool evalFunction_2_lt(const std::string & name, Value & arg0, Value & arg1, Value & result, std::false_type);
bool evalFunction_2_lt(const std::string & name, Value & arg0, Value & arg1, Value & result, std::true_type);
void evalNumericRange(const std::string & str, Value & 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 & 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 split(const std::string & str, const char delimeter);
template static bool isNumber(const std::string & str, T * num = 0);
template static T stringToNumber(const std::string & str);
template 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 ParserXd;
typedef Parser ParserXf;
typedef Parser ParserXi;
//----------------------------------------
// Function definitions.
//----------------------------------------
template
Parser::Parser() :
mOperators1("+-*/^()[]="),
mOperators2(".+.-.*./.^"),
mCacheChunkedExpressions(false)
{
// Coefficient-wise operations.
mFunctions.push_back("abs");
mFunctions.push_back("sqrt");
mFunctions.push_back("exp");
mFunctions.push_back("log");
mFunctions.push_back("log10");
mFunctions.push_back("sin");
mFunctions.push_back("cos");
mFunctions.push_back("tan");
mFunctions.push_back("asin");
mFunctions.push_back("acos");
// Matrix reduction operations.
mFunctions.push_back("trace");
mFunctions.push_back("norm");
mFunctions.push_back("size");
if (has_operator_lt::value) {
mFunctions.push_back("min");
mFunctions.push_back("max");
mFunctions.push_back("absmax");
}
mFunctions.push_back("mean");
mFunctions.push_back("sum");
mFunctions.push_back("prod");
// Matrix operations.
mFunctions.push_back("transpose");
mFunctions.push_back("conjugate");
mFunctions.push_back("adjoint");
// Matrix initializers.
mFunctions.push_back("zeros");
mFunctions.push_back("ones");
mFunctions.push_back("eye");
}
template
Value Parser::eval(const std::string & expression)
{
#ifdef DEBUG
# ifdef EIGENLAB_DEBUG
std::cout 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 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 field == "^" || op->field == ".^")) {
if(lhs->type == VARIABLE) {
if(!isVariable(lhs->field))
throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + lhs->field + "'.");
lhs->value.setShared(mVariables[lhs->field]);
}
if(rhs->type == VARIABLE) {
if(!isVariable(rhs->field))
throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + rhs->field + "'.");
rhs->value.setShared(mVariables[rhs->field]);
}
if(rhs->value.matrix().size() == 1) {
lhs->value.local() = lhs->value.matrix().array().pow(rhs->value.matrix()(0, 0));
lhs->value.mapLocal();
lhs->type = VALUE;
} else if(lhs->value.matrix().size() == 1) {
typename Derived::Scalar temp = lhs->value.matrix()(0, 0);
lhs->value.local().resize(rhs->value.matrix().rows(), rhs->value.matrix().cols());
for(size_t row = 0; row < size_t(rhs->value.matrix().rows()); row++) {
for(size_t col = 0; col < size_t(rhs->value.matrix().cols()); col++)
lhs->value.local()(row, col) = pow(temp, rhs->value.matrix()(row, col));
}
lhs->value.mapLocal();
lhs->type = VALUE;
} else if(op->field == ".^" && lhs->value.matrix().rows() == rhs->value.matrix().rows() && lhs->value.matrix().cols() == rhs->value.matrix().cols()) {
lhs->value.local().resize(rhs->value.matrix().rows(), rhs->value.matrix().cols());
for(size_t row = 0; row < size_t(rhs->value.matrix().rows()); row++) {
for(size_t col = 0; col < size_t(rhs->value.matrix().cols()); col++)
lhs->value.local()(row, col) = pow(lhs->value.matrix()(row, col), rhs->value.matrix()(row, col));
}
lhs->value.mapLocal();
lhs->type = VALUE;
} else {
throw std::runtime_error("Invalid operand dimensions for operation '" + lhs->field + op->field + rhs->field + "'.");
}
chunks.erase(op, rhs + 1);
op = lhs + 1;
rhs = (op != chunks.end()) ? op + 1 : op;
#ifdef DEBUG
# ifdef EIGENLAB_DEBUG
operationPerformed = true;
# endif
#endif
} else {
lhs = op;
op = rhs;
rhs++;
}
}
#ifdef DEBUG
# ifdef EIGENLAB_DEBUG
if(operationPerformed) { std::cout field == "*" || op->field == "/" || op->field == ".*" || op->field == "./")) {
if(lhs->type == VARIABLE) {
if(!isVariable(lhs->field))
throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + lhs->field + "'.");
lhs->value.setShared(mVariables[lhs->field]);
}
if(rhs->type == VARIABLE) {
if(!isVariable(rhs->field))
throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + rhs->field + "'.");
rhs->value.setShared(mVariables[rhs->field]);
}
if(rhs->value.matrix().size() == 1) {
if(lhs->value.isLocal()) {
if(op->field == "*" || op->field == ".*")
lhs->value.local().array() *= rhs->value.matrix()(0, 0);
else // if(op->field == "/" || op->field == "./")
lhs->value.local().array() /= rhs->value.matrix()(0, 0);
} else {
if(op->field == "*" || op->field == ".*")
lhs->value.local() = lhs->value.matrix().array() * rhs->value.matrix()(0, 0);
else // if(op->field == "/" || op->field == "./")
lhs->value.local() = lhs->value.matrix().array() / rhs->value.matrix()(0, 0);
lhs->value.mapLocal();
lhs->type = VALUE;
}
} else if(lhs->value.matrix().size() == 1) {
typename Derived::Scalar temp = lhs->value.matrix()(0, 0);
if(op->field == "*" || op->field == ".*")
lhs->value.local() = rhs->value.matrix().array() * temp;
else // if(op->field == "/" || op->field == "./")
lhs->value.local() = Derived::Constant(rhs->value.matrix().rows(), rhs->value.matrix().cols(), temp).array() / rhs->value.matrix().array();
lhs->value.mapLocal();
lhs->type = VALUE;
} else if((op->field == ".*" || op->field == "./") && lhs->value.matrix().rows() == rhs->value.matrix().rows() && lhs->value.matrix().cols() == rhs->value.matrix().cols()) {
if(lhs->value.isLocal()) {
if(op->field == ".*")
lhs->value.local().array() *= rhs->value.matrix().array();
else // if(op->field == "./")
lhs->value.local().array() /= rhs->value.matrix().array();
} else {
if(op->field == ".*")
lhs->value.local() = lhs->value.matrix().array() * rhs->value.matrix().array();
else // if(op->field == "./")
lhs->value.local() = lhs->value.matrix().array() / rhs->value.matrix().array();
lhs->value.mapLocal();
lhs->type = VALUE;
}
} else if(op->field == "*" && lhs->value.matrix().cols() == rhs->value.matrix().rows()) {
if(lhs->value.isLocal()) {
lhs->value.local() *= rhs->value.matrix();
lhs->value.mapLocal();
} else {
lhs->value.local() = lhs->value.matrix() * rhs->value.matrix();
lhs->value.mapLocal();
lhs->type = VALUE;
}
} else {
throw std::runtime_error("Invalid operand dimensions for operation '" + lhs->field + op->field + rhs->field + "'.");
}
chunks.erase(op, rhs + 1);
op = lhs + 1;
rhs = (op != chunks.end()) ? op + 1 : op;
#ifdef DEBUG
# ifdef EIGENLAB_DEBUG
operationPerformed = true;
# endif
#endif
} else {
lhs = op;
op = rhs;
rhs++;
}
}
#ifdef DEBUG
# ifdef EIGENLAB_DEBUG
if(operationPerformed) { std::cout field == "+" || op->field == "-" || op->field == ".+" || op->field == ".-")) {
if(lhs->type == VARIABLE) {
if(!isVariable(lhs->field))
throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + lhs->field + "'.");
lhs->value.setShared(mVariables[lhs->field]);
}
if(rhs->type == VARIABLE) {
if(!isVariable(rhs->field))
throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + rhs->field + "'.");
rhs->value.setShared(mVariables[rhs->field]);
}
if(rhs->value.matrix().size() == 1) {
if(lhs->value.isLocal()) {
if(op->field == "+" || op->field == ".+")
lhs->value.local().array() += rhs->value.matrix()(0, 0);
else // if(op->field == "-" || op->field == ".-")
lhs->value.local().array() -= rhs->value.matrix()(0, 0);
} else {
if(op->field == "+" || op->field == ".+")
lhs->value.local() = lhs->value.matrix().array() + rhs->value.matrix()(0, 0);
else // if(op->field == "-" || op->field == ".-")
lhs->value.local() = lhs->value.matrix().array() - rhs->value.matrix()(0, 0);
lhs->value.mapLocal();
lhs->type = VALUE;
}
} else if(lhs->value.matrix().size() == 1) {
typename Derived::Scalar temp = lhs->value.matrix()(0, 0);
if(op->field == "+" || op->field == ".+")
lhs->value.local() = rhs->value.matrix().array() + temp;
else // if(op->field == "-" || op->field == ".-")
lhs->value.local() = Derived::Constant(rhs->value.matrix().rows(), rhs->value.matrix().cols(), temp).array() - rhs->value.matrix().array();
lhs->value.mapLocal();
lhs->type = VALUE;
} else if(lhs->value.matrix().rows() == rhs->value.matrix().rows() && lhs->value.matrix().cols() == rhs->value.matrix().cols()) {
if(lhs->value.isLocal()) {
if(op->field == "+" || op->field == ".+")
lhs->value.local().array() += rhs->value.matrix().array();
else // if(op->field == "-" || op->field == ".-")
lhs->value.local().array() -= rhs->value.matrix().array();
} else {
if(op->field == "+" || op->field == ".+")
lhs->value.local() = lhs->value.matrix().array() + rhs->value.matrix().array();
else // if(op->field == "-" || op->field == ".-")
lhs->value.local() = lhs->value.matrix().array() - rhs->value.matrix().array();
lhs->value.mapLocal();
lhs->type = VALUE;
}
} else {
throw std::runtime_error("Invalid operand dimensions for operation '" + lhs->field + op->field + rhs->field + "'.");
}
chunks.erase(op, rhs + 1);
op = lhs + 1;
rhs = (op != chunks.end()) ? op + 1 : op;
#ifdef DEBUG
# ifdef EIGENLAB_DEBUG
operationPerformed = true;
# endif
#endif
} else {
lhs = op;
op = rhs;
rhs++;
}
}
#ifdef DEBUG
# ifdef EIGENLAB_DEBUG
if(operationPerformed) { std::cout 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