EigenLab.h
Go to the documentation of this file.
1 // --*- Mode: C++; c-basic-offset:4; indent-tabs-mode:t; tab-width:4 -*--
2 // EigenLab
3 // Version: 1.0.0
4 // Author: Dr. Marcel Paz Goldschen-Ohm
5 // Email: marcel.goldschen@gmail.com
6 // Copyright (c) 2015 by Dr. Marcel Paz Goldschen-Ohm.
7 // Licence: MIT
8 //----------------------------------------
9 
10 #ifndef EigenLab_H
11 #define EigenLab_H
12 
13 #include <Eigen/Dense>
14 #include <iomanip>
15 #include <type_traits>
16 #include <map>
17 #include <sstream>
18 #include <stdexcept>
19 #include <string>
20 #include <vector>
21 
22 // Define both DEBUG and EIGENLAB_DEBUG for step-by-step equation parsing printouts.
23 #ifndef DEBUG
24 //# define DEBUG
25 #endif
26 
27 #ifndef EIGENLAB_DEBUG
28 //# define EIGENLAB_DEBUG
29 #endif
30 
31 #ifdef DEBUG
32 # include <iostream>
33 #endif
34 
35 namespace EigenLab
36 {
37  //----------------------------------------
38  // A wrapper for a matrix whose data is either stored locally or shared.
39  //
40  // Template typename Derived can be any dynamically sized matrix type supported by Eigen.
41  //
42  // !!! matrix() promises to ALWAYS return a map to the matrix data whether it's
43  // stored locally or externally in some shared memory.
44  //
45  // !!! local() provides direct access to the local data, but this data is
46  // ONLY valid when isLocal() is true. In most cases, you're best off
47  // accessing the matrix data via matrix() instead.
48  //----------------------------------------
49  template <typename Derived = Eigen::MatrixXd>
50  class Value
51  {
52  private:
53  // Local matrix data.
54  Derived mLocal;
55 
56  // Map to shared matrix data (map points to mLocal if the data is local).
57  // !!! This map promises to ALWAYS point to the matrix data whether it's
58  // stored locally in mLocal or externally in some shared memory.
59  Eigen::Map<Derived> mShared;
60 
61  // Flag indicating whether the local data is being used.
62  bool mIsLocal;
63 
64  public:
65  // Access mapped data (whether its local or shared).
66  // !!! matrix() promises to ALWAYS return a map to the matrix data whether it's
67  // stored locally in mLocal or externally in some shared memory.
68  inline Eigen::Map<Derived> & matrix() { return mShared; }
69  inline const Eigen::Map<Derived> & matrix() const { return mShared; }
70 
71  // Access local data.
72  // !!! WARNING! This data is ONLY valid if isLocal() is true.
73  // !!! WARNING! If you change the local data via this method, you MUST call mapLocal() immediately afterwards.
74  // In most cases, you're best off accessing the matrix data via matrix() instead.
75  inline Derived & local() { return mLocal; }
76  inline const Derived & local() const { return mLocal; }
77 
78  // Is mapped data local?
79  inline bool isLocal() const { return mIsLocal; }
80 
81  // Set mapped data to point to local data.
82  inline void mapLocal() { new (& mShared) Eigen::Map<Derived>(mLocal.data(), mLocal.rows(), mLocal.cols()); mIsLocal = true; }
83 
84  // Copy shared data to local data (if needed).
85  inline void copySharedToLocal() { if(!isLocal()) { mLocal = mShared; mapLocal(); } }
86 
87  // Set local data.
88  Value() : mLocal(1, 1), mShared(mLocal.data(), mLocal.rows(), mLocal.cols()), mIsLocal(true) {}
89  Value(const typename Derived::Scalar s) : mLocal(Derived::Constant(1, 1, s)), mShared(mLocal.data(), mLocal.rows(), mLocal.cols()), mIsLocal(true) {}
90  Value(const Derived & mat) : mLocal(mat), mShared(mLocal.data(), mLocal.rows(), mLocal.cols()), mIsLocal(true) {}
91  inline void setLocal(const typename Derived::Scalar s) { mLocal = Derived::Constant(1, 1, s); mapLocal(); }
92  inline void setLocal(const Eigen::MatrixBase<Derived> & mat) { mLocal = mat; mapLocal(); }
93  inline void setLocal(const Value & val) { mLocal = val.matrix(); mapLocal(); }
94  inline void setLocal(const typename Derived::Scalar * data, size_t rows = 1, size_t cols = 1) { setShared(data, rows, cols); copySharedToLocal(); }
95  inline Value & operator = (const typename Derived::Scalar s) { setLocal(s); return (* this); }
96  inline Value & operator = (const Derived & mat) { setLocal(mat); return (* this); }
97 
98  // Set shared data.
99  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) {}
100  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; }
101  inline void setShared(const Derived & mat) { setShared(mat.data(), mat.rows(), mat.cols()); }
102  inline void setShared(const Value & val) { setShared(val.matrix().data(), val.matrix().rows(), val.matrix().cols()); }
103 
104  // Set to local or shared data dependent on whether val maps its own local data or some other shared data.
105  Value(const Value & val) : mLocal(1, 1), mShared(mLocal.data(), mLocal.rows(), mLocal.cols()) { (* this) = val; }
106  inline Value & operator = (const Value & val) { if(val.isLocal()) { setLocal(val); } else { setShared(val); } return (* this); }
107  };
111 
112  // check if a class has a comparison operator (ie. std::complex does not)
113  template<typename T>
115  {
116  template<class U>
117  static auto test(U*) -> decltype(std::declval<U>() < std::declval<U>());
118  template<typename>
119  static auto test(...) -> std::false_type;
120  using type = typename std::is_same<bool, decltype(test<T>(0))>::type;
121  };
122  template<typename T>
123  struct has_operator_lt : has_operator_lt_impl<T>::type {};
124 
125  //----------------------------------------
126  // Equation parser.
127  //
128  // Template typename Derived can be any dynamically sized matrix type supported by Eigen.
129  //----------------------------------------
130  template <typename Derived = Eigen::MatrixXd>
131  class Parser
132  {
133  public:
134  // A map to hold named values.
135  typedef std::map<std::string, Value<Derived> > ValueMap;
136 
137  private:
138  // Variables are stored in a map of named values.
139  ValueMap mVariables;
140 
141  // Operator symbols and function names used by the parser.
142  std::string mOperators1, mOperators2;
143  std::vector<std::string> mFunctions;
144 
145  // Expressions are parsed by first splitting them into chunks.
146  struct Chunk {
147  std::string field;
148  int type;
150  int row0, col0, rows, cols;
151  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) {}
152  };
153  enum ChunkType { VALUE = 0, VARIABLE, OPERATOR, FUNCTION };
154  typedef std::vector<Chunk> ChunkArray;
155  typedef typename Derived::Index Index;
157  std::map<std::string, ChunkArray> mCachedChunkedExpressions;
158 
159  public:
160  // Access to named variables.
161  // !!! WARNING! var(name) will create the variable name if it does not already exist.
162  inline ValueMap & vars() { return mVariables; }
163  inline Value<Derived> & var(const std::string & name) { return mVariables[name]; }
164 
165  // Check if a variable exists.
166  inline bool hasVar(const std::string & name) { return isVariable(name); }
167 
168  // Delete a variable.
169  inline void clearVar(const std::string & name) { typename ValueMap::iterator it = mVariables.find(name); if(it != mVariables.end()) mVariables.erase(it); }
170 
171  // Expression chunk caching.
172  inline bool cacheExpressions() const { return mCacheChunkedExpressions; }
173  inline void setCacheExpressions(bool b) { mCacheChunkedExpressions = b; }
174  inline void clearCachedExpressions() { mCachedChunkedExpressions.clear(); }
175 
176  Parser();
177  ~Parser() { clearCachedExpressions(); }
178 
179  // Evaluate an expression and return the result in a value wrapper.
180  Value<Derived> eval(const std::string & expression);
181 
182  private:
183  void splitEquationIntoChunks(const std::string & expression, ChunkArray & chunks, std::string & code);
184  std::string::const_iterator findClosingBracket(const std::string & str, const std::string::const_iterator openingBracket, const char closingBracket) const;
185  std::vector<std::string> splitArguments(const std::string & str, const char delimeter) const;
186  void evalIndexRange(const std::string & str, int * first, int * last, int numIndices);
187  void evalMatrixExpression(const std::string & str, Value<Derived> & mat);
188  void evalFunction(const std::string & name, std::vector<std::string> & args, Value<Derived> & result);
189  bool evalFunction_1_lt(const std::string & name, Value<Derived> & arg, Value<Derived> & result, std::false_type);
190  bool evalFunction_1_lt(const std::string & name, Value<Derived> & arg, Value<Derived> & result, std::true_type);
191  bool evalFunction_2_lt(const std::string & name, Value<Derived> & arg0, Value<Derived> & arg1, Value<Derived> & result, std::false_type);
192  bool evalFunction_2_lt(const std::string & name, Value<Derived> & arg0, Value<Derived> & arg1, Value<Derived> & result, std::true_type);
193 
194  void evalNumericRange(const std::string & str, Value<Derived> & mat);
195  inline bool isVariable(const std::string & name) const { return mVariables.count(name) > 0; }
196  inline bool isOperator(const char c) const { return (std::find(mOperators1.begin(), mOperators1.end(), c) != mOperators1.end()); }
197  bool isOperator(const std::string & str) const;
198  inline bool isFunction(const std::string & str) const { return (std::find(mFunctions.begin(), mFunctions.end(), str) != mFunctions.end()); }
199  void evalIndices(ChunkArray & chunks);
200  void evalNegations(ChunkArray & chunks);
201  void evalPowers(ChunkArray & chunks);
202  void evalMultiplication(ChunkArray & chunks);
203  void evalAddition(ChunkArray & chunks);
204  void evalAssignment(ChunkArray & chunks);
205 #ifdef DEBUG
206 # ifdef EIGENLAB_DEBUG
207  void printChunks(ChunkArray & chunks, size_t maxRows = 2, size_t maxCols = 2, int precision = 0);
208  void printVars(size_t maxRows = 2, size_t maxCols = 2, int precision = 0);
209  std::string textRepresentation(Value<Derived> & val, size_t maxRows = 2, size_t maxCols = 2, int precision = 0);
210 # endif
211 #endif
212 
213  public:
214  static std::string trim(const std::string & str);
215  static std::vector<std::string> split(const std::string & str, const char delimeter);
216  template <typename T> static bool isNumber(const std::string & str, T * num = 0);
217  template <typename T> static T stringToNumber(const std::string & str);
218  template <typename T> static std::string numberToString(T num, int precision = 0);
219 #ifdef DEBUG
220  void test_w_lt(size_t & numFails, typename Derived::Scalar & s, Derived & a34, Derived & b34, Derived & c43, Derived & v, std::true_type);
221  void test_w_lt(size_t & numFails, typename Derived::Scalar & s, Derived & a34, Derived & b34, Derived & c43, Derived & v, std::false_type);
222  size_t test();
223 #endif
224  };
228 
229  //----------------------------------------
230  // Function definitions.
231  //----------------------------------------
232  template <typename Derived>
234  mOperators1("+-*/^()[]="),
235  mOperators2(".+.-.*./.^"),
236  mCacheChunkedExpressions(false)
237  {
238  // Coefficient-wise operations.
239  mFunctions.push_back("abs");
240  mFunctions.push_back("sqrt");
241  mFunctions.push_back("square");
242  mFunctions.push_back("exp");
243  mFunctions.push_back("log");
244  mFunctions.push_back("log10");
245  mFunctions.push_back("sin");
246  mFunctions.push_back("cos");
247  mFunctions.push_back("tan");
248  mFunctions.push_back("asin");
249  mFunctions.push_back("acos");
250 
251  // Matrix reduction operations.
252  mFunctions.push_back("trace");
253  mFunctions.push_back("norm");
254  mFunctions.push_back("size");
256  mFunctions.push_back("min");
257  mFunctions.push_back("minOfFinites");
258  mFunctions.push_back("max");
259  mFunctions.push_back("maxOfFinites");
260  mFunctions.push_back("absmax");
261  mFunctions.push_back("cwiseMin");
262  mFunctions.push_back("cwiseMax");
263  }
264  mFunctions.push_back("mean");
265  mFunctions.push_back("meanOfFinites");
266  mFunctions.push_back("sum");
267  mFunctions.push_back("sumOfFinites");
268  mFunctions.push_back("prod");
269  mFunctions.push_back("numberOfFinites");
270 
271  // Matrix operations.
272  mFunctions.push_back("transpose");
273  mFunctions.push_back("conjugate");
274  mFunctions.push_back("adjoint");
275 
276  // Matrix initializers.
277  mFunctions.push_back("zeros");
278  mFunctions.push_back("ones");
279  mFunctions.push_back("eye");
280  }
281 
282  template <typename Derived>
283  Value<Derived> Parser<Derived>::eval(const std::string & expression)
284  {
285 #ifdef DEBUG
286 # ifdef EIGENLAB_DEBUG
287  std::cout << "---" << std::endl;
288  std::cout << "EXPRESSION: " << expression << std::endl;
289 # endif
290 #endif
291  ChunkArray chunks;
292  std::string code;
293  splitEquationIntoChunks(trim(expression), chunks, code);
294  evalIndices(chunks);
295  evalNegations(chunks);
296  evalPowers(chunks);
297  evalMultiplication(chunks);
298  evalAddition(chunks);
299  evalAssignment(chunks);
300  if(chunks.size() != 1)
301  throw std::runtime_error("Failed to reduce expression '" + expression + "' to a single value.");
302 #ifdef DEBUG
303 # ifdef EIGENLAB_DEBUG
304  std::cout << "---" << std::endl;
305 # endif
306 #endif
307  if(chunks[0].type == VARIABLE) {
308  if(!isVariable(chunks[0].field))
309  throw std::runtime_error("Unknown variable '" + chunks[0].field + "'.");
310  return mVariables[chunks[0].field];
311  }
312  return chunks[0].value;
313  }
314 
315  template <typename Derived>
316  void Parser<Derived>::splitEquationIntoChunks(const std::string & expression, ChunkArray & chunks, std::string & code)
317  {
318  if(cacheExpressions()) {
319  if(mCachedChunkedExpressions.count(expression) > 0) {
320  chunks = mCachedChunkedExpressions.at(expression);
321 #ifdef DEBUG
322 # ifdef EIGENLAB_DEBUG
323  std::cout << "CACHED CHUNKS: "; printChunks(chunks); std::cout << std::endl;
324 # endif
325 #endif
326  return;
327  }
328  }
329 
330  for(std::string::const_iterator it = expression.begin(); it != expression.end();)
331  {
332  int prevType = (chunks.size() ? chunks.back().type : -1);
333  char ci = (* it);
334  if(isspace(ci)) {
335  // Ignore whitespace.
336  it++;
337  } else if(ci == '(' && (prevType == VALUE || prevType == VARIABLE)) {
338  // Index group.
339  std::string::const_iterator jt = findClosingBracket(expression, it, ')');
340  if(jt == expression.end())
341  throw std::runtime_error("Missing closing bracket for '" + std::string(it, jt) + "'.");
342  std::string field = std::string(it + 1, jt); // Outer brackets stripped.
343  if(prevType == VARIABLE) {
344  if(!isVariable(chunks.back().field))
345  throw std::runtime_error("Unknown variable '" + chunks.back().field + "'.");
346  chunks.back().value.setShared(var(chunks.back().field));
347  }
348  int first, last;
349  int rows = int(chunks.back().value.matrix().rows());
350  int cols = int(chunks.back().value.matrix().cols());
351  std::vector<std::string> args = splitArguments(field, ',');
352  if(args.size() == 1) {
353  if(cols == 1) {
354  evalIndexRange(args[0], & first, & last, rows);
355  chunks.back().row0 = first;
356  chunks.back().col0 = 0;
357  chunks.back().rows = last + 1 - first;//(last == -1 ? int(chunks.back().value.matrix().rows()) : last + 1) - first;
358  chunks.back().cols = 1;
359  } else if(rows == 1) {
360  evalIndexRange(args[0], & first, & last, cols);
361  chunks.back().row0 = 0;
362  chunks.back().col0 = first;
363  chunks.back().rows = 1;
364  chunks.back().cols = last + 1 - first;//(last == -1 ? int(chunks.back().value.matrix().cols()) : last + 1) - first;
365  } else {
366  throw std::runtime_error("Missing row or column indices for '(" + chunks.back().field + "(" + field + ")'.");
367  }
368  } else if(args.size() == 2) {
369  evalIndexRange(args[0], & first, & last, rows);
370  chunks.back().row0 = first;
371  chunks.back().rows = last + 1 - first;//(last == -1 ? int(chunks.back().value.matrix().rows()) : last + 1) - first;
372  evalIndexRange(args[1], & first, & last, cols);
373  chunks.back().col0 = first;
374  chunks.back().cols = last + 1 - first;//(last == -1 ? int(chunks.back().value.matrix().cols()) : last + 1) - first;
375  } else {
376  throw std::runtime_error("Invalid index expression '" + chunks.back().field + "(" + field + ")'.");
377  }
378  code += "i";
379  it = jt + 1;
380  } else if(ci == '(' && prevType == FUNCTION) {
381  // Function argument group.
382  std::string::const_iterator jt = findClosingBracket(expression, it, ')');
383  if(jt == expression.end())
384  throw std::runtime_error("Missing closing bracket for '" + std::string(it, jt) + "'.");
385  std::string field = std::string(it + 1, jt); // Outer brackets stripped.
386  std::vector<std::string> args = splitArguments(field, ',');
387  evalFunction(chunks.back().field, args, chunks.back().value);
388  chunks.back().field += "(" + field + ")";
389  chunks.back().type = VALUE;
390  code += (chunks.back().value.matrix().size() == 1 ? "n" : "M");
391  it = jt + 1;
392  } else if(ci == '(') {
393  // Recursively evaluate group to a single value.
394  std::string::const_iterator jt = findClosingBracket(expression, it, ')');
395  if(jt == expression.end())
396  throw std::runtime_error("Missing closing bracket for '" + std::string(it, jt) + "'.");
397  std::string field = std::string(it + 1, jt); // Outer brackets stripped.
398  chunks.push_back(Chunk(field, prevType = VALUE, eval(field)));
399  code += (chunks.back().value.matrix().size() == 1 ? "n" : "M");
400  it = jt + 1;
401  } else if(ci == '[') {
402  // Evaluate matrix.
403  if(prevType == VALUE || prevType == VARIABLE)
404  throw std::runtime_error("Invalid operation '" + chunks.back().field + std::string(1, ci) + "'.");
405  std::string::const_iterator jt = findClosingBracket(expression, it, ']');
406  if(jt == expression.end())
407  throw std::runtime_error("Missing closing bracket for '" + std::string(it, jt) + "'.");
408  std::string field = std::string(it + 1, jt); // Outer brackets stripped.
409  chunks.push_back(Chunk("[" + field + "]", prevType = VALUE));
410  evalMatrixExpression(field, chunks.back().value);
411  code += (chunks.back().value.matrix().size() == 1 ? "n" : "M");
412  it = jt + 1;
413  } else if(it + 1 != expression.end() && isOperator(std::string(it, it + 2))) {
414  // Double character operator.
415  std::string field = std::string(it, it + 2);
416  chunks.push_back(Chunk(field, prevType = OPERATOR));
417  code += field;
418  it += 2;
419  } else if(isOperator(ci)) {
420  // Single character operator.
421  std::string field = std::string(1, ci);
422  chunks.push_back(Chunk(field, prevType = OPERATOR));
423  code += field;
424  it++;
425  } else {
426  // Non-operator: value range, number, function, or variable name.
427  std::string::const_iterator jt = it + 1;
428  // accept fp-strings, ie [+-]
429  unsigned char state = 1;
430  for(std::string::const_iterator kt = it; state && kt != expression.end(); kt++) {
431  unsigned char token;
432  if (*kt == ' ') token = 0;
433  else if (*kt == '+' || *kt == '-') token = 1;
434  else if (isdigit(*kt)) token = 2;
435  else if (*kt == '.') token = 3;
436  else if (*kt == 'e' || *kt == 'E') token = 4;
437  else break;
438  const static char nextstate[9][5] = {{0},
439  {1, 2, 3, 4, 0},
440  {0, 0, 3, 4, 0},
441  {0, 0, 3, 5, 6},
442  {0, 0, 5, 0, 0},
443  {0, 0, 5, 0, 6},
444  {0, 7, 8, 0, 0},
445  {0, 0, 8, 0, 0},
446  {0, 0, 8, 0, 0}};
447  //WARN("state=" << (int)state << " token(" << *kt << ")=" << (int)token
448  // << " nextstate = " << (int)nextstate[state][token] << "\n");
449  state = nextstate[state][token];
450  if (state == 8) jt = kt;
451  }
452  for(; jt != expression.end(); jt++) {
453  if(isOperator(*jt) || (jt + 1 != expression.end() && isOperator(std::string(jt, jt + 2))))
454  break;
455  }
456  std::string field = trim(std::string(it, jt));
457  if(prevType == VALUE || prevType == VARIABLE)
458  throw std::runtime_error("Invalid operation '" + chunks.back().field + field + "'.");
459  double num;
460  if(field.find(":") != std::string::npos) {
461  // Numeric range.
462  chunks.push_back(Chunk(field, prevType = VALUE));
463  evalNumericRange(field, chunks.back().value);
464  code += (chunks.back().value.matrix().size() == 1 ? "n" : "M");
465  } else if(isNumber<double>(field, & num)) {
466  // Number.
467  chunks.push_back(Chunk(field, prevType = VALUE, Value<Derived>(num)));
468  code += "n";
469  } else if(isVariable(field)) {
470  // Local variable.
471  chunks.push_back(Chunk(field, prevType = VARIABLE));
472  code += (mVariables[field].matrix().size() == 1 ? "vn" : "vM");
473  } else if(isFunction(field)) {
474  // Function.
475  chunks.push_back(Chunk(field, prevType = FUNCTION));
476  } else {
477  // New undefined variable.
478  chunks.push_back(Chunk(field, prevType = VARIABLE));
479  code += "vU";
480  }
481  it = jt;
482  }
483  } // it
484 #ifdef DEBUG
485 # ifdef EIGENLAB_DEBUG
486  std::cout << "CHUNKS: "; printChunks(chunks); std::cout << std::endl;
487  std::cout << "CODE: " << code << std::endl;
488 # endif
489 #endif
490  if(cacheExpressions())
491  mCachedChunkedExpressions[expression] = chunks;
492  }
493 
494  template <typename Derived>
495  std::string::const_iterator Parser<Derived>::findClosingBracket(const std::string & str, const std::string::const_iterator openingBracket, const char closingBracket) const
496  {
497  int depth = 1;
498  std::string::const_iterator it = openingBracket + 1;
499  for(; it != str.end(); it++) {
500  if((*it) == (*openingBracket)) depth++;
501  else if((*it) == closingBracket) depth--;
502  if(depth == 0) return it;
503  }
504  return str.end();
505  }
506 
507  template <typename Derived>
508  std::vector<std::string> Parser<Derived>::splitArguments(const std::string & str, const char delimeter) const
509  {
510  std::vector<std::string> args;
511  std::string::const_iterator i0 = str.begin();
512  for(std::string::const_iterator it = str.begin(); it != str.end(); it++) {
513  if((*it) == '(') it = findClosingBracket(str, it, ')');
514  else if((*it) == '[') it = findClosingBracket(str, it, ']');
515  else if((*it) == delimeter) {
516  args.push_back(trim(std::string(i0, it)));
517  i0 = it + 1;
518  }
519  }
520  args.push_back(std::string(i0, str.end()));
521  return args;
522  }
523 
524  template <typename Derived>
525  void Parser<Derived>::evalIndexRange(const std::string & str, int * first, int * last, int numIndices)
526  {
527  if(str.empty())
528  throw std::runtime_error("Empty index range.");
529  ValueXi valuei;
530  ParserXi parseri;
531  size_t pos;
532  for(std::string::const_iterator it = str.begin(); it != str.end(); it++) {
533  if((*it) == ':') {
534  std::string firstStr = trim(std::string(str.begin(), it));
535  std::string lastStr = trim(std::string(it + 1, str.end()));
536  if(firstStr.empty() && lastStr.empty()) {
537  (* first) = 0;
538  (* last) = numIndices - 1;
539  return;
540  }
541  if(firstStr.empty() || lastStr.empty())
542  throw std::runtime_error("Missing indices for '" + str + "'.");
543 
544  pos = firstStr.find("end");
545  if(pos != std::string::npos) {
546  firstStr = firstStr.substr(0, pos) + numberToString<int>(numIndices - 1) + firstStr.substr(pos + 3);
547  }
548  pos = lastStr.find("end");
549  if(pos != std::string::npos) {
550  lastStr = lastStr.substr(0, pos) + numberToString<int>(numIndices - 1) + lastStr.substr(pos + 3);
551  }
552 
553  valuei = parseri.eval(firstStr);
554  if(valuei.matrix().size() != 1)
555  throw std::runtime_error("Invalid indices '" + str + "'.");
556  (* first) = valuei.matrix()(0, 0);
557 
558  valuei = parseri.eval(lastStr);
559  if(valuei.matrix().size() != 1)
560  throw std::runtime_error("Invalid indices '" + str + "'.");
561  (* last) = valuei.matrix()(0, 0);
562 
563  return;
564  }
565  }
566  std::string firstStr = str;
567 
568  pos = firstStr.find("end");
569  if(pos != std::string::npos) {
570  firstStr = firstStr.substr(0, pos) + numberToString<int>(numIndices - 1) + firstStr.substr(pos + 3);
571  }
572 
573  valuei = parseri.eval(firstStr);
574  if(valuei.matrix().size() != 1)
575  throw std::runtime_error("Invalid index '" + str + "'.");
576  (* first) = valuei.matrix()(0, 0);
577  (* last) = (* first);
578  }
579 
580  template <typename Derived>
581  void Parser<Derived>::evalMatrixExpression(const std::string & str, Value<Derived> & mat)
582  {
583  // !!! Expression may NOT include outer brackets, although brackets for individual rows are OK.
584  std::vector<std::string> rows = splitArguments(str, ';');
585  std::vector<std::vector<typename Derived::Scalar> > temp;
586  Value<Derived> submatrix;
587  size_t row0 = 0, col0 = 0, nrows = 0, ncols = 0;
588  size_t pos;
589  for(size_t i = 0; i < rows.size(); i++) {
590  // Strip row brackets if they exist.
591  if(rows[i][0] == '[' && rows[i].back() == ']') rows[i] = rows[i].substr(1, int(rows[i].size()) - 2);
592  std::vector<std::string> cols = splitArguments(rows[i], ',');
593  col0 = 0;
594  ncols = 0;
595  for(size_t j = 0; j < cols.size(); j++) {
596  pos = cols[j].find(":");
597  if(pos != std::string::npos) {
598  std::string firstStr = cols[j].substr(0, pos);
599  std::string lastStr = cols[j].substr(pos + 1);
600  pos = lastStr.find(":");
601  if(pos != std::string::npos) {
602  std::string stepStr = lastStr.substr(0, pos);
603  lastStr = lastStr.substr(pos + 1);
604  if(lastStr.find(":") != std::string::npos)
605  throw std::runtime_error("Invalid matrix definition '[" + str + "]'. Invalid range '" + cols[j] + "'.");
606  // first:step:last
607  Value<Derived> first = eval(firstStr);
608  Value<Derived> step = eval(stepStr);
609  Value<Derived> last = eval(lastStr);
610  if(first.matrix().size() != 1 || step.matrix().size() != 1 || last.matrix().size() != 1)
611  throw std::runtime_error("Invalid matrix definition '[" + str + "]'. Invalid range '" + cols[j] + "'.");
612  typename Derived::RealScalar sfirst = std::real(first.matrix()(0));
613  typename Derived::RealScalar sstep = std::real(step.matrix()(0));
614  typename Derived::RealScalar slast = std::real(last.matrix()(0));
615  if(sfirst == slast) {
616  submatrix.local().setConstant(1, 1, sfirst);
617  submatrix.mapLocal();
618  } else if((slast - sfirst >= 0 && sstep > 0) || (slast - sfirst <= 0 && sstep < 0)) {
619  int n = floor((slast - sfirst) / sstep) + 1;
620  submatrix.local().resize(1, n);
621  for(int k = 0; k < n; ++k)
622  submatrix.local()(0, k) = sfirst + k * sstep;
623  submatrix.mapLocal();
624  } else {
625  throw std::runtime_error("Invalid matrix definition '[" + str + "]'. Invalid range '" + cols[j] + "'.");
626  }
627  } else {
628  // first:last => first:1:last
629  Value<Derived> first = eval(firstStr);
630  Value<Derived> last = eval(lastStr);
631  if(first.matrix().size() != 1 || last.matrix().size() != 1)
632  throw std::runtime_error("Invalid matrix definition '[" + str + "]'. Invalid range '" + cols[j] + "'.");
633  typename Derived::RealScalar sfirst = std::real(first.matrix()(0));
634  typename Derived::RealScalar slast = std::real(last.matrix()(0));
635  if(sfirst == slast) {
636  submatrix.local().setConstant(1, 1, sfirst);
637  submatrix.mapLocal();
638  } else if(slast - sfirst >= 0) {
639  int n = floor(slast - sfirst) + 1;
640  submatrix.local().resize(1, n);
641  for(int k = 0; k < n; ++k)
642  submatrix.local()(0, k) = sfirst + k;
643  submatrix.mapLocal();
644  } else {
645  throw std::runtime_error("Invalid matrix definition '[" + str + "]'. Invalid range '" + cols[j] + "'.");
646  }
647  }
648  } else {
649  submatrix = eval(cols[j]);
650  }
651  if(j > 0 && size_t(submatrix.matrix().cols()) != nrows)
652  throw std::runtime_error("Invalid matrix definition '[" + str + "]'. Successive column entries '" + cols[int(j) - 1] + "' and '" + cols[j] + "' do not have the same number of rows.");
653  nrows = submatrix.matrix().rows();
654  ncols += submatrix.matrix().cols();
655  temp.resize(row0 + submatrix.matrix().rows());
656  for(size_t row = 0; row < size_t(submatrix.matrix().rows()); row++) {
657  temp[row0 + row].resize(col0 + submatrix.matrix().cols());
658  for(size_t col = 0; col < size_t(submatrix.matrix().cols()); col++)
659  temp[row0 + row][col0 + col] = submatrix.matrix()(row, col);
660  }
661  col0 += submatrix.matrix().cols();
662  }
663  if(row0 > 0 && ncols != temp[int(row0) - 1].size())
664  throw std::runtime_error("Invalid matrix definition '[" + str + "]'. Successive row entries '" + rows[int(i) - 1] + "' and '" + rows[i] + "' do not have the same number of columns.");
665  row0 += nrows;
666  }
667  nrows = temp.size();
668  if(nrows == 0) return;
669  ncols = temp[0].size();
670  mat.setLocal(Derived(nrows, ncols));
671  for(size_t row = 0; row < nrows; row++) {
672  for(size_t col = 0; col < ncols; col++)
673  mat.matrix()(row, col) = temp[row][col];
674  }
675  mat.mapLocal();
676  }
677 
678  template <typename Derived>
679  bool Parser<Derived>::evalFunction_1_lt(const std::string & name, Value<Derived> & arg, Value<Derived> & result, std::true_type)
680  {
681  if(name == "min") {
682  result.setLocal(arg.matrix().minCoeff());
683  return true;
684  } else if(name == "minOfFinites") {
685  result.setLocal(arg.matrix().minCoeffOfFinites());
686  return true;
687  } else if(name == "max") {
688  result.setLocal(arg.matrix().maxCoeff());
689  return true;
690  } else if(name == "maxOfFinites") {
691  result.setLocal(arg.matrix().maxCoeffOfFinites());
692  return true;
693  } else if(name == "absmax") {
694  typename Derived::Scalar minimum = arg.matrix().minCoeff();
695  typename Derived::Scalar maximum = arg.matrix().maxCoeff();
696  result.setLocal(std::abs(maximum) >= std::abs(minimum) ? maximum : minimum);
697  return true;
698  }
699  return false;
700  }
701 
702  template <typename Derived>
703  bool Parser<Derived>::evalFunction_1_lt(const std::string &/*name*/, Value<Derived> &/*arg0*/, Value<Derived> &/*result*/, std::false_type)
704  {
705  return false;
706  }
707 
708  template <typename Derived>
709  bool Parser<Derived>::evalFunction_2_lt(const std::string & name, Value<Derived> & arg0, Value<Derived> & arg1, Value<Derived> & result, std::true_type)
710  {
711  if(name == "min") {
712  if(arg1.matrix().size() != 1)
713  throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
714  int dim = floor(std::real(arg1.matrix()(0, 0)));
715  if((dim != 0 && dim != 1) || dim != std::real(arg1.matrix()(0, 0)))
716  throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
717  if(dim == 0) {
718  result.local() = arg0.matrix().colwise().minCoeff();
719  result.mapLocal();
720  return true;
721  } else if(dim == 1) {
722  result.local() = arg0.matrix().rowwise().minCoeff();
723  result.mapLocal();
724  return true;
725  }
726  } else if(name == "max") {
727  if(arg1.matrix().size() != 1)
728  throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
729  int dim = floor(std::real(arg1.matrix()(0, 0)));
730  if((dim != 0 && dim != 1) || dim != std::real(arg1.matrix()(0, 0)))
731  throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
732  if(dim == 0) {
733  result.local() = arg0.matrix().colwise().maxCoeff();
734  result.mapLocal();
735  return true;
736  } else if(dim == 1) {
737  result.local() = arg0.matrix().rowwise().maxCoeff();
738  result.mapLocal();
739  return true;
740  }
741  } else if(name == "absmax") {
742  if(arg1.matrix().size() != 1)
743  throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
744  int dim = floor(std::real(arg1.matrix()(0, 0)));
745  if((dim != 0 && dim != 1) || dim != std::real(arg1.matrix()(0, 0)))
746  throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
747  if(dim == 0) {
748  result.local() = arg0.matrix().colwise().maxCoeff();
749  result.mapLocal();
750  Derived minimum = arg0.matrix().colwise().minCoeff();
751  for(size_t i = 0; i < size_t(result.matrix().size()); i++) {
752  if(std::abs(result.matrix()(i)) < std::abs(minimum(i)))
753  result.matrix()(i) = minimum(i);
754  }
755  return true;
756  } else if(dim == 1) {
757  result.local() = arg0.matrix().rowwise().maxCoeff();
758  result.mapLocal();
759  Derived minimum = arg0.matrix().rowwise().minCoeff();
760  for(size_t i = 0; i < size_t(result.matrix().size()); i++) {
761  if(std::abs(result.matrix()(i)) < std::abs(minimum(i)))
762  result.matrix()(i) = minimum(i);
763  }
764  return true;
765  }
766  } else if (name == "cwiseMin") {
767  if (arg1.matrix().size() == 1) {
768  typename Derived::Scalar arg1scalar = arg1.matrix()(0, 0);
769  Derived arg1matrix = Derived::Constant(arg0.matrix().rows(), arg0.matrix().cols(), arg1scalar);
770  result.local() = arg0.matrix().cwiseMin(arg1matrix);
771  result.mapLocal();
772  return true;
773  } else if (arg0.matrix().cols() == arg1.matrix().cols() && arg0.matrix().rows() == arg1.matrix().rows()) {
774  result.local() = arg0.matrix().cwiseMin(arg1.matrix());
775  result.mapLocal();
776  return true;
777  } else {
778  throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
779  }
780  } else if (name == "cwiseMax") {
781  if (arg1.matrix().size() == 1) {
782  typename Derived::Scalar arg1scalar = arg1.matrix()(0, 0);
783  Derived arg1matrix = Derived::Constant(arg0.matrix().rows(), arg0.matrix().cols(), arg1scalar);
784  result.local() = arg0.matrix().cwiseMax(arg1matrix);
785  result.mapLocal();
786  return true;
787  } else if (arg0.matrix().cols() == arg1.matrix().cols() && arg0.matrix().rows() == arg1.matrix().rows()) {
788  result.local() = arg0.matrix().cwiseMax(arg1.matrix());
789  result.mapLocal();
790  return true;
791  } else {
792  throw std::runtime_error("Invalid dimension argument for function '" + name + "(...)'.");
793  }
794  }
795  return false;
796  }
797 
798  template <typename Derived>
799  bool Parser<Derived>::evalFunction_2_lt(const std::string &/*name*/, Value<Derived> &/*arg0*/, Value<Derived> &/*arg1*/, Value<Derived> &/*result*/, std::false_type)
800  {
801  return false;
802  }
803 
804  template <typename Derived>
805  void Parser<Derived>::evalFunction(const std::string & name, std::vector<std::string> & args, Value<Derived> & result)
806  {
807  if(args.size() == 1) {
808  Value<Derived> arg = eval(args[0]);
809  if(name == "abs") {
810  result.local() = arg.matrix().array().abs().template cast<typename Derived::Scalar>();
811  result.mapLocal();
812  return;
813  } else if(name == "sqrt") {
814  result.local() = arg.matrix().array().sqrt();
815  result.mapLocal();
816  return;
817  } else if(name == "square") {
818  result.local() = arg.matrix().array().square();
819  result.mapLocal();
820  return;
821  } else if(name == "exp") {
822  result.local() = arg.matrix().array().exp();
823  result.mapLocal();
824  return;
825  } else if(name == "log") {
826  result.local() = arg.matrix().array().log();
827  result.mapLocal();
828  return;
829  } else if(name == "log10") {
830  result.local() = arg.matrix().array().log();
831  result.local() *= (1.0 / log(10));
832  result.mapLocal();
833  return;
834  } else if(name == "sin") {
835  result.local() = arg.matrix().array().sin();
836  result.mapLocal();
837  return;
838  } else if(name == "cos") {
839  result.local() = arg.matrix().array().cos();
840  result.mapLocal();
841  return;
842  } else if(name == "tan") {
843  result.local() = arg.matrix().array().tan();
844  result.mapLocal();
845  return;
846  } else if(name == "acos") {
847  result.local() = arg.matrix().array().acos();
848  result.mapLocal();
849  return;
850  } else if(name == "asin") {
851  result.local() = arg.matrix().array().asin();
852  result.mapLocal();
853  return;
854  } else if(name == "trace") {
855  result.setLocal(arg.matrix().trace());
856  return;
857  } else if(name == "norm") {
858  result.setLocal(arg.matrix().norm());
859  return;
860  } else if (evalFunction_1_lt(name, arg, result, has_operator_lt<typename Derived::Scalar>())) {
861  return;
862  } else if(name == "mean") {
863  result.setLocal(arg.matrix().mean());
864  return;
865  } else if(name == "meanOfFinites") {
866  result.setLocal(arg.matrix().meanOfFinites());
867  return;
868  } else if(name == "sum") {
869  result.setLocal(arg.matrix().sum());
870  return;
871  } else if(name == "sumOfFinites") {
872  result.setLocal(arg.matrix().sumOfFinites());
873  return;
874  } else if(name == "prod") {
875  result.setLocal(arg.matrix().prod());
876  return;
877  } else if(name == "numberOfFinites") {
878  result.setLocal(arg.matrix().numberOfFinites());
879  return;
880  } else if(name == "transpose") {
881  result.local() = arg.matrix().transpose();
882  result.mapLocal();
883  return;
884  } else if(name == "conjugate") {
885  result.local() = arg.matrix().conjugate();
886  result.mapLocal();
887  return;
888  } else if(name == "adjoint") {
889  result.local() = arg.matrix().adjoint();
890  result.mapLocal();
891  return;
892  } else if(name == "zeros") {
893  if(arg.matrix().size() != 1)
894  throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + ")'.");
895  int N = floor(std::real(arg.matrix()(0, 0)));
896  if(N <= 0 || N != std::real(arg.matrix()(0, 0)))
897  throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + ")'.");
898  result.local() = Derived::Zero(N, N);
899  result.mapLocal();
900  return;
901  } else if(name == "ones") {
902  if(arg.matrix().size() != 1)
903  throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + ")'.");
904  int N = floor(std::real(arg.matrix()(0, 0)));
905  if(N <= 0 || N != std::real(arg.matrix()(0, 0)))
906  throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + ")'.");
907  result.local() = Derived::Ones(N, N);
908  result.mapLocal();
909  return;
910  } else if(name == "eye") {
911  if(arg.matrix().size() != 1)
912  throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + ")'.");
913  int N = floor(std::real(arg.matrix()(0, 0)));
914  if(N <= 0 || N != std::real(arg.matrix()(0, 0)))
915  throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + ")'.");
916  result.local() = Derived::Identity(N, N);
917  result.mapLocal();
918  return;
919  }
920  else {
921  throw std::runtime_error("Invalid function '" + name + "(" + args[0] + ")'.");
922  }
923  } else if(args.size() == 2) {
924  Value<Derived> arg0 = eval(args[0]);
925  Value<Derived> arg1 = eval(args[1]);
926  if(name == "size") {
927  if(arg1.matrix().size() != 1)
928  throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
929  int dim = floor(std::real(arg1.matrix()(0, 0)));
930  if((dim != 0 && dim != 1) || dim != std::real(arg1.matrix()(0, 0)))
931  throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
932  if(dim == 0) {
933  result.setLocal((typename Derived::Scalar) arg0.matrix().rows());
934  return;
935  } else if(dim == 1) {
936  result.setLocal((typename Derived::Scalar) arg0.matrix().cols());
937  return;
938  }
939  } else if (evalFunction_2_lt(name, arg0, arg1, result, has_operator_lt<typename Derived::Scalar>())) {
940  return;
941  } else if(name == "mean") {
942  if(arg1.matrix().size() != 1)
943  throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
944  int dim = floor(std::real(arg1.matrix()(0, 0)));
945  if((dim != 0 && dim != 1) || dim != std::real(arg1.matrix()(0, 0)))
946  throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
947  if(dim == 0) {
948  result.local() = arg0.matrix().colwise().mean();
949  result.mapLocal();
950  return;
951  } else if(dim == 1) {
952  result.local() = arg0.matrix().rowwise().mean();
953  result.mapLocal();
954  return;
955  }
956  } else if(name == "sum") {
957  if(arg1.matrix().size() != 1)
958  throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
959  int dim = floor(std::real(arg1.matrix()(0, 0)));
960  if((dim != 0 && dim != 1) || dim != std::real(arg1.matrix()(0, 0)))
961  throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
962  if(dim == 0) {
963  result.local() = arg0.matrix().colwise().sum();
964  result.mapLocal();
965  return;
966  } else if(dim == 1) {
967  result.local() = arg0.matrix().rowwise().sum();
968  result.mapLocal();
969  return;
970  }
971  } else if(name == "prod") {
972  if(arg1.matrix().size() != 1)
973  throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
974  int dim = floor(std::real(arg1.matrix()(0, 0)));
975  if((dim != 0 && dim != 1) || dim != std::real(arg1.matrix()(0, 0)))
976  throw std::runtime_error("Invalid dimension argument for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
977  if(dim == 0) {
978  result.local() = arg0.matrix().colwise().prod();
979  result.mapLocal();
980  return;
981  } else if(dim == 1) {
982  result.local() = arg0.matrix().rowwise().prod();
983  result.mapLocal();
984  return;
985  }
986  } else if(name == "zeros") {
987  if((arg0.matrix().size() != 1) || (arg1.matrix().size() != 1))
988  throw std::runtime_error("Invalid dimension arguments for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
989  int rows = floor(std::real(arg0.matrix()(0, 0)));
990  int cols = floor(std::real(arg1.matrix()(0, 0)));
991  if(rows <= 0 || cols <= 0 || rows != std::real(arg0.matrix()(0, 0)) || cols != std::real(arg1.matrix()(0, 0)))
992  throw std::runtime_error("Invalid dimension arguments for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
993  result.local() = Derived::Zero(rows, cols);
994  result.mapLocal();
995  return;
996  } else if(name == "ones") {
997  if((arg0.matrix().size() != 1) || (arg1.matrix().size() != 1))
998  throw std::runtime_error("Invalid dimension arguments for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
999  int rows = floor(std::real(arg0.matrix()(0, 0)));
1000  int cols = floor(std::real(arg1.matrix()(0, 0)));
1001  if(rows <= 0 || cols <= 0 || rows != std::real(arg0.matrix()(0, 0)) || cols != std::real(arg1.matrix()(0, 0)))
1002  throw std::runtime_error("Invalid dimension arguments for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
1003  result.local() = Derived::Ones(rows, cols);
1004  result.mapLocal();
1005  return;
1006  } else if(name == "eye") {
1007  if((arg0.matrix().size() != 1) || (arg1.matrix().size() != 1))
1008  throw std::runtime_error("Invalid dimension arguments for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
1009  int rows = floor(std::real(arg0.matrix()(0, 0)));
1010  int cols = floor(std::real(arg1.matrix()(0, 0)));
1011  if(rows <= 0 || cols <= 0 || rows != std::real(arg0.matrix()(0, 0)) || cols != std::real(arg1.matrix()(0, 0)))
1012  throw std::runtime_error("Invalid dimension arguments for function '" + name + "(" + args[0] + "," + args[1] + ")'.");
1013  result.local() = Derived::Identity(rows, cols);
1014  result.mapLocal();
1015  return;
1016  } else {
1017  throw std::runtime_error("Invalid function '" + name + "(" + args[0] + "," + args[1] + ")'.");
1018  }
1019  }
1020  std::string argsStr = "(";
1021  for(size_t i = 0; i < args.size(); i++) {
1022  if(i > 0) argsStr += ",";
1023  argsStr += args[i];
1024  }
1025  argsStr += ")";
1026  throw std::runtime_error("Invalid function/arguments for '" + name + argsStr + "'.");
1027  }
1028 
1029  template <typename Derived>
1030  void Parser<Derived>::evalNumericRange(const std::string & str, Value<Derived> & mat)
1031  {
1032  size_t pos = str.find(":");
1033  if(pos == std::string::npos)
1034  throw std::runtime_error("Invalid numeric range '" + str + "'.");
1035  size_t pos2 = str.substr(pos + 1).find(":");
1036  if(pos2 == std::string::npos) {
1037  // first:last
1038  std::string firstStr = str.substr(0, pos);
1039  std::string lastStr = str.substr(pos + 1);
1040  Value<Derived> first = eval(firstStr);
1041  Value<Derived> last = eval(lastStr);
1042  if(first.matrix().size() != 1 || last.matrix().size() != 1)
1043  throw std::runtime_error("Invalid numeric range '" + str + "'.");
1044  typename Derived::RealScalar sfirst = std::real(first.matrix()(0,0));
1045  typename Derived::RealScalar slast = std::real(last.matrix()(0,0));
1046  if(sfirst > slast)
1047  throw std::runtime_error("Invalid numeric range '" + str + "'. Must not reverse.");
1048  int n = 1 + floor(slast - sfirst);
1049  mat.local().resize(1, n);
1050  for(int i = 0; i < n; i++)
1051  mat.local()(0, i) = sfirst + i;
1052  mat.mapLocal();
1053  } else {
1054  // first:step:last
1055  pos2 += pos + 1;
1056  std::string firstStr = str.substr(0, pos);
1057  std::string stepStr = str.substr(pos + 1, pos2 - pos - 1);
1058  std::string lastStr = str.substr(pos2 + 1);
1059  Value<Derived> first = eval(firstStr);
1060  Value<Derived> step = eval(stepStr);
1061  Value<Derived> last = eval(lastStr);
1062  if(first.matrix().size() != 1 || step.matrix().size() != 1 || last.matrix().size() != 1)
1063  throw std::runtime_error("Invalid numeric range '" + str + "'.");
1064  typename Derived::RealScalar sfirst = std::real(first.matrix()(0, 0));
1065  typename Derived::RealScalar sstep = std::real(step.matrix()(0, 0));
1066  typename Derived::RealScalar slast = std::real(last.matrix()(0, 0));
1067  if(sfirst == slast) {
1068  mat = sfirst;
1069  } else if(sfirst < slast && sstep > 0) {
1070  int n = 1 + floor((slast - sfirst) / sstep);
1071  mat.local().resize(1, n);
1072  for(int i = 0; i < n; i++)
1073  mat.local()(0, i) = sfirst + i * sstep;
1074  mat.mapLocal();
1075  } else if(sfirst > slast && sstep < 0) {
1076  int n = 1 + floor((slast - sfirst) / sstep);
1077  mat.local().resize(1, n);
1078  for(int i = 0; i < n; i++)
1079  mat.local()(0, i) = sfirst + i * sstep;
1080  mat.mapLocal();
1081  } else {
1082  throw std::runtime_error("Invalid numeric range '" + str + "'.");
1083  }
1084  }
1085  }
1086 
1087  template <typename Derived>
1088  bool Parser<Derived>::isOperator(const std::string & str) const
1089  {
1090  if(str.size() == 1)
1091  return isOperator(str[0]);
1092  else if(str.size() == 2) {
1093  size_t pos = mOperators2.find(str);
1094  return (pos != std::string::npos && pos % 2 == 0);
1095  }
1096  return false;
1097  }
1098 
1099  template <typename Derived>
1101  {
1102 #ifdef DEBUG
1103 # ifdef EIGENLAB_DEBUG
1104  bool operationPerformed = false;
1105 # endif
1106 #endif
1107  for(typename ChunkArray::iterator it = chunks.begin(); it != chunks.end(); it++) {
1108  if(it->row0 != -1 && (it->type == VALUE || (it->type == VARIABLE && (it + 1 == chunks.end() || (it + 1)->type != OPERATOR || (it + 1)->field != "=")))) {
1109  if(it->type == VALUE) {
1110  Derived temp = it->value.local().block(it->row0, it->col0, it->rows, it->cols);
1111  it->value.local() = temp;
1112  it->value.mapLocal();
1113  } else { //if(it->type == VARIABLE) {
1114  if(!isVariable(it->field))
1115  throw std::runtime_error("Attempted indexing into uninitialized variable '" + it->field + "'.");
1116  it->value.local() = mVariables[it->field].matrix().block(it->row0, it->col0, it->rows, it->cols);
1117  it->value.mapLocal();
1118  it->type = VALUE;
1119  }
1120  it->row0 = -1;
1121  it->col0 = -1;
1122  it->rows = -1;
1123  it->cols = -1;
1124 #ifdef DEBUG
1125 # ifdef EIGENLAB_DEBUG
1126  operationPerformed = true;
1127 # endif
1128 #endif
1129  }
1130  }
1131 #ifdef DEBUG
1132 # ifdef EIGENLAB_DEBUG
1133  if(operationPerformed) { std::cout << "i: "; printChunks(chunks); std::cout << std::endl; }
1134 # endif
1135 #endif
1136  }
1137 
1138  template <typename Derived>
1140  {
1141 #ifdef DEBUG
1142 # ifdef EIGENLAB_DEBUG
1143  bool operationPerformed = false;
1144 # endif
1145 #endif
1146  if(chunks.size() < 2) return;
1147  typename ChunkArray::iterator lhs = chunks.begin(), op = chunks.begin(), rhs = op + 1;
1148  for(; lhs != chunks.end() && op != chunks.end() && rhs != chunks.end();) {
1149  if(op->type == OPERATOR && op->field == "-" && (op == chunks.begin() || (lhs->type != VALUE && lhs->type != VARIABLE)) && (rhs->type == VALUE || rhs->type == VARIABLE)) {
1150  if(rhs->type == VALUE)
1151  rhs->value.matrix().array() *= -1;
1152  else if(rhs->type == VARIABLE) {
1153  if(!isVariable(rhs->field))
1154  throw std::runtime_error("Attempted operation '" + op->field + rhs->field + "' on uninitialized variable '" + rhs->field + "'.");
1155  rhs->value.local() = mVariables[rhs->field].matrix().array() * -1;
1156  rhs->value.mapLocal();
1157  rhs->type = VALUE;
1158  }
1159  lhs = chunks.erase(op);
1160  op = (lhs != chunks.end()) ? lhs + 1 : lhs;
1161  rhs = (op != chunks.end()) ? op + 1 : op;
1162 #ifdef DEBUG
1163 # ifdef EIGENLAB_DEBUG
1164  operationPerformed = true;
1165 # endif
1166 #endif
1167  } else {
1168  lhs = op;
1169  op = rhs;
1170  rhs++;
1171  }
1172  }
1173 #ifdef DEBUG
1174 # ifdef EIGENLAB_DEBUG
1175  if(operationPerformed) { std::cout << "-: "; printChunks(chunks); std::cout << std::endl; }
1176 # endif
1177 #endif
1178  }
1179 
1180  template <typename Derived>
1182  {
1183 #ifdef DEBUG
1184 # ifdef EIGENLAB_DEBUG
1185  bool operationPerformed = false;
1186 # endif
1187 #endif
1188  if(chunks.size() < 3) return;
1189  typename ChunkArray::iterator lhs = chunks.begin(), op = lhs + 1, rhs = op + 1;
1190  for(; lhs != chunks.end() && op != chunks.end() && rhs != chunks.end();) {
1191  if((op->type == OPERATOR) && (op->field == "^" || op->field == ".^")) {
1192  if(lhs->type == VARIABLE) {
1193  if(!isVariable(lhs->field))
1194  throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + lhs->field + "'.");
1195  lhs->value.setShared(mVariables[lhs->field]);
1196  }
1197  if(rhs->type == VARIABLE) {
1198  if(!isVariable(rhs->field))
1199  throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + rhs->field + "'.");
1200  rhs->value.setShared(mVariables[rhs->field]);
1201  }
1202  if(rhs->value.matrix().size() == 1) {
1203  lhs->value.local() = lhs->value.matrix().array().pow(rhs->value.matrix()(0, 0));
1204  lhs->value.mapLocal();
1205  lhs->type = VALUE;
1206  } else if(lhs->value.matrix().size() == 1) {
1207  typename Derived::Scalar temp = lhs->value.matrix()(0, 0);
1208  lhs->value.local().resize(rhs->value.matrix().rows(), rhs->value.matrix().cols());
1209  for(size_t row = 0; row < size_t(rhs->value.matrix().rows()); row++) {
1210  for(size_t col = 0; col < size_t(rhs->value.matrix().cols()); col++)
1211  lhs->value.local()(row, col) = pow(temp, rhs->value.matrix()(row, col));
1212  }
1213  lhs->value.mapLocal();
1214  lhs->type = VALUE;
1215  } else if(op->field == ".^" && lhs->value.matrix().rows() == rhs->value.matrix().rows() && lhs->value.matrix().cols() == rhs->value.matrix().cols()) {
1216  lhs->value.local().resize(rhs->value.matrix().rows(), rhs->value.matrix().cols());
1217  for(size_t row = 0; row < size_t(rhs->value.matrix().rows()); row++) {
1218  for(size_t col = 0; col < size_t(rhs->value.matrix().cols()); col++)
1219  lhs->value.local()(row, col) = pow(lhs->value.matrix()(row, col), rhs->value.matrix()(row, col));
1220  }
1221  lhs->value.mapLocal();
1222  lhs->type = VALUE;
1223  } else {
1224  throw std::runtime_error("Invalid operand dimensions for operation '" + lhs->field + op->field + rhs->field + "'.");
1225  }
1226  chunks.erase(op, rhs + 1);
1227  op = lhs + 1;
1228  rhs = (op != chunks.end()) ? op + 1 : op;
1229 #ifdef DEBUG
1230 # ifdef EIGENLAB_DEBUG
1231  operationPerformed = true;
1232 # endif
1233 #endif
1234  } else {
1235  lhs = op;
1236  op = rhs;
1237  rhs++;
1238  }
1239  }
1240 #ifdef DEBUG
1241 # ifdef EIGENLAB_DEBUG
1242  if(operationPerformed) { std::cout << "^: "; printChunks(chunks); std::cout << std::endl; }
1243 # endif
1244 #endif
1245  }
1246 
1247  template <typename Derived>
1249  {
1250 #ifdef DEBUG
1251 # ifdef EIGENLAB_DEBUG
1252  bool operationPerformed = false;
1253 # endif
1254 #endif
1255  if(chunks.size() < 3) return;
1256  typename ChunkArray::iterator lhs = chunks.begin(), op = lhs + 1, rhs = op + 1;
1257  for(; lhs != chunks.end() && op != chunks.end() && rhs != chunks.end();) {
1258  if((op->type == OPERATOR) && (op->field == "*" || op->field == "/" || op->field == ".*" || op->field == "./")) {
1259  if(lhs->type == VARIABLE) {
1260  if(!isVariable(lhs->field))
1261  throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + lhs->field + "'.");
1262  lhs->value.setShared(mVariables[lhs->field]);
1263  }
1264  if(rhs->type == VARIABLE) {
1265  if(!isVariable(rhs->field))
1266  throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + rhs->field + "'.");
1267  rhs->value.setShared(mVariables[rhs->field]);
1268  }
1269  if(rhs->value.matrix().size() == 1) {
1270  if(lhs->value.isLocal()) {
1271  if(op->field == "*" || op->field == ".*")
1272  lhs->value.local().array() *= rhs->value.matrix()(0, 0);
1273  else // if(op->field == "/" || op->field == "./")
1274  lhs->value.local().array() /= rhs->value.matrix()(0, 0);
1275  } else {
1276  if(op->field == "*" || op->field == ".*")
1277  lhs->value.local() = lhs->value.matrix().array() * rhs->value.matrix()(0, 0);
1278  else // if(op->field == "/" || op->field == "./")
1279  lhs->value.local() = lhs->value.matrix().array() / rhs->value.matrix()(0, 0);
1280  lhs->value.mapLocal();
1281  lhs->type = VALUE;
1282  }
1283  } else if(lhs->value.matrix().size() == 1) {
1284  typename Derived::Scalar temp = lhs->value.matrix()(0, 0);
1285  if(op->field == "*" || op->field == ".*")
1286  lhs->value.local() = rhs->value.matrix().array() * temp;
1287  else // if(op->field == "/" || op->field == "./")
1288  lhs->value.local() = Derived::Constant(rhs->value.matrix().rows(), rhs->value.matrix().cols(), temp).array() / rhs->value.matrix().array();
1289  lhs->value.mapLocal();
1290  lhs->type = VALUE;
1291  } else if((op->field == ".*" || op->field == "./") && lhs->value.matrix().rows() == rhs->value.matrix().rows() && lhs->value.matrix().cols() == rhs->value.matrix().cols()) {
1292  if(lhs->value.isLocal()) {
1293  if(op->field == ".*")
1294  lhs->value.local().array() *= rhs->value.matrix().array();
1295  else // if(op->field == "./")
1296  lhs->value.local().array() /= rhs->value.matrix().array();
1297  } else {
1298  if(op->field == ".*")
1299  lhs->value.local() = lhs->value.matrix().array() * rhs->value.matrix().array();
1300  else // if(op->field == "./")
1301  lhs->value.local() = lhs->value.matrix().array() / rhs->value.matrix().array();
1302  lhs->value.mapLocal();
1303  lhs->type = VALUE;
1304  }
1305  } else if(op->field == "*" && lhs->value.matrix().cols() == rhs->value.matrix().rows()) {
1306  if(lhs->value.isLocal()) {
1307  lhs->value.local() *= rhs->value.matrix();
1308  lhs->value.mapLocal();
1309  } else {
1310  lhs->value.local() = lhs->value.matrix() * rhs->value.matrix();
1311  lhs->value.mapLocal();
1312  lhs->type = VALUE;
1313  }
1314  } else {
1315  throw std::runtime_error("Invalid operand dimensions for operation '" + lhs->field + op->field + rhs->field + "'.");
1316  }
1317  chunks.erase(op, rhs + 1);
1318  op = lhs + 1;
1319  rhs = (op != chunks.end()) ? op + 1 : op;
1320 #ifdef DEBUG
1321 # ifdef EIGENLAB_DEBUG
1322  operationPerformed = true;
1323 # endif
1324 #endif
1325  } else {
1326  lhs = op;
1327  op = rhs;
1328  rhs++;
1329  }
1330  }
1331 #ifdef DEBUG
1332 # ifdef EIGENLAB_DEBUG
1333  if(operationPerformed) { std::cout << "*: "; printChunks(chunks); std::cout << std::endl; }
1334 # endif
1335 #endif
1336  }
1337 
1338  template <typename Derived>
1340  {
1341 #ifdef DEBUG
1342 # ifdef EIGENLAB_DEBUG
1343  bool operationPerformed = false;
1344 # endif
1345 #endif
1346  if(chunks.size() < 3) return;
1347  typename ChunkArray::iterator lhs = chunks.begin(), op = lhs + 1, rhs = op + 1;
1348  for(; lhs != chunks.end() && op != chunks.end() && rhs != chunks.end();) {
1349  if((op->type == OPERATOR) && (op->field == "+" || op->field == "-" || op->field == ".+" || op->field == ".-")) {
1350  if(lhs->type == VARIABLE) {
1351  if(!isVariable(lhs->field))
1352  throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + lhs->field + "'.");
1353  lhs->value.setShared(mVariables[lhs->field]);
1354  }
1355  if(rhs->type == VARIABLE) {
1356  if(!isVariable(rhs->field))
1357  throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + rhs->field + "'.");
1358  rhs->value.setShared(mVariables[rhs->field]);
1359  }
1360  if(rhs->value.matrix().size() == 1) {
1361  if(lhs->value.isLocal()) {
1362  if(op->field == "+" || op->field == ".+")
1363  lhs->value.local().array() += rhs->value.matrix()(0, 0);
1364  else // if(op->field == "-" || op->field == ".-")
1365  lhs->value.local().array() -= rhs->value.matrix()(0, 0);
1366  } else {
1367  if(op->field == "+" || op->field == ".+")
1368  lhs->value.local() = lhs->value.matrix().array() + rhs->value.matrix()(0, 0);
1369  else // if(op->field == "-" || op->field == ".-")
1370  lhs->value.local() = lhs->value.matrix().array() - rhs->value.matrix()(0, 0);
1371  lhs->value.mapLocal();
1372  lhs->type = VALUE;
1373  }
1374  } else if(lhs->value.matrix().size() == 1) {
1375  typename Derived::Scalar temp = lhs->value.matrix()(0, 0);
1376  if(op->field == "+" || op->field == ".+")
1377  lhs->value.local() = rhs->value.matrix().array() + temp;
1378  else // if(op->field == "-" || op->field == ".-")
1379  lhs->value.local() = Derived::Constant(rhs->value.matrix().rows(), rhs->value.matrix().cols(), temp).array() - rhs->value.matrix().array();
1380  lhs->value.mapLocal();
1381  lhs->type = VALUE;
1382  } else if(lhs->value.matrix().rows() == rhs->value.matrix().rows() && lhs->value.matrix().cols() == rhs->value.matrix().cols()) {
1383  if(lhs->value.isLocal()) {
1384  if(op->field == "+" || op->field == ".+")
1385  lhs->value.local().array() += rhs->value.matrix().array();
1386  else // if(op->field == "-" || op->field == ".-")
1387  lhs->value.local().array() -= rhs->value.matrix().array();
1388  } else {
1389  if(op->field == "+" || op->field == ".+")
1390  lhs->value.local() = lhs->value.matrix().array() + rhs->value.matrix().array();
1391  else // if(op->field == "-" || op->field == ".-")
1392  lhs->value.local() = lhs->value.matrix().array() - rhs->value.matrix().array();
1393  lhs->value.mapLocal();
1394  lhs->type = VALUE;
1395  }
1396  } else {
1397  throw std::runtime_error("Invalid operand dimensions for operation '" + lhs->field + op->field + rhs->field + "'.");
1398  }
1399  chunks.erase(op, rhs + 1);
1400  op = lhs + 1;
1401  rhs = (op != chunks.end()) ? op + 1 : op;
1402 #ifdef DEBUG
1403 # ifdef EIGENLAB_DEBUG
1404  operationPerformed = true;
1405 # endif
1406 #endif
1407  } else {
1408  lhs = op;
1409  op = rhs;
1410  rhs++;
1411  }
1412  }
1413 #ifdef DEBUG
1414 # ifdef EIGENLAB_DEBUG
1415  if(operationPerformed) { std::cout << "+: "; printChunks(chunks); std::cout << std::endl; }
1416 # endif
1417 #endif
1418  }
1419 
1420  template <typename Derived>
1422  {
1423 #ifdef DEBUG
1424 # ifdef EIGENLAB_DEBUG
1425  bool operationPerformed = false;
1426 # endif
1427 #endif
1428  if(chunks.size() < 3) return;
1429  typename ChunkArray::iterator rhs = chunks.end() - 1, op = rhs - 1, lhs = op - 1;
1430  for(; op != chunks.begin() && rhs != chunks.begin();) {
1431  if(op->type == OPERATOR && op->field == "=" && (lhs->type == VALUE || lhs->type == VARIABLE) && (rhs->type == VALUE || rhs->type == VARIABLE)) {
1432  if(rhs->type == VARIABLE) {
1433  if(!isVariable(rhs->field))
1434  throw std::runtime_error("Attempted operation '" + lhs->field + op->field + rhs->field + "' on uninitialized variable '" + rhs->field + "'.");
1435  rhs->value.setShared(mVariables[rhs->field]);
1436  }
1437  if(lhs->type == VALUE) {
1438  lhs->value.local() = rhs->value.matrix();
1439  lhs->value.mapLocal();
1440  } else { //if(lhs->type == VARIABLE) {
1441  if(isVariable(lhs->field)) {
1442  lhs->value.setShared(mVariables[lhs->field]);
1443  if(lhs->row0 == -1) {
1444  if(lhs->value.matrix().rows() == rhs->value.matrix().rows() && lhs->value.matrix().cols() == rhs->value.matrix().cols()) {
1445  lhs->value.matrix() = rhs->value.matrix();
1446  } else {
1447  mVariables[lhs->field].local() = rhs->value.matrix();
1448  mVariables[lhs->field].mapLocal();
1449  }
1450  } else { //if(lhs->row0 != -1) {
1451  lhs->value.matrix().block(lhs->row0, lhs->col0, lhs->rows, lhs->cols) = rhs->value.matrix();
1452  }
1453  } else {
1454  mVariables[lhs->field].local() = rhs->value.matrix();
1455  mVariables[lhs->field].mapLocal();
1456  }
1457  }
1458  rhs = chunks.erase(op, rhs + 1);
1459  op = (rhs != chunks.begin()) ? rhs - 1 : rhs;
1460  if (op != chunks.begin()) lhs = op - 1;
1461 #ifdef DEBUG
1462 # ifdef EIGENLAB_DEBUG
1463  operationPerformed = true;
1464 # endif
1465 #endif
1466  } else {
1467  rhs = op;
1468  op = lhs;
1469  lhs--;
1470  }
1471  }
1472 #ifdef DEBUG
1473 # ifdef EIGENLAB_DEBUG
1474  if(operationPerformed) { std::cout << "=: "; printChunks(chunks); std::cout << std::endl; }
1475 # endif
1476 #endif
1477  }
1478 
1479 #ifdef DEBUG
1480 # ifdef EIGENLAB_DEBUG
1481  template <typename Derived>
1482  void Parser<Derived>::printChunks(ChunkArray & chunks, size_t maxRows, size_t maxCols, int precision)
1483  {
1484  std::cout << "__";
1485  for(typename ChunkArray::iterator it = chunks.begin(); it != chunks.end(); it++) {
1486  switch(it->type) {
1487  case VALUE:
1488  std::cout << textRepresentation(it->value, maxRows, maxCols, precision);
1489  if(it->row0 != -1)
1490  std::cout << "(" << it->row0 << ":" << it->row0 + it->rows - 1 << "," << it->col0 << ":" << it->col0 + it->cols - 1 << ")";
1491  break;
1492  case VARIABLE:
1493  std::cout << it->field;// << "=" << textRepresentation(mVariables[it->field], maxRows, maxCols, precision);
1494  if(it->row0 != -1)
1495  std::cout << "(" << it->row0 << ":" << it->row0 + it->rows - 1 << "," << it->col0 << ":" << it->col0 + it->cols - 1 << ")";
1496  break;
1497  case OPERATOR:
1498  std::cout << it->field;
1499  break;
1500  case FUNCTION:
1501  std::cout << "f()=" << it->field;
1502  break;
1503  }
1504  std::cout << "__";
1505  }
1506  }
1507 
1508  template <typename Derived>
1509  void Parser<Derived>::printVars(size_t maxRows, size_t maxCols, int precision)
1510  {
1511  for(typename ValueMap::iterator it = mVariables.begin(); it != mVariables.end(); it++)
1512  std::cout << it->first << " (" << it->second.matrix().rows() << "x" << it->second.matrix().cols() << ") = " << textRepresentation(it->second, maxRows, maxCols, precision) << std::endl;
1513  }
1514 
1515  template <typename Derived>
1516  std::string Parser<Derived>::textRepresentation(Value<Derived> & val, size_t maxRows, size_t maxCols, int precision)
1517  {
1518  if(val.matrix().size() == 1)
1519  return numberToString<typename Derived::Scalar>(val.matrix()(0, 0), precision);
1520  else {
1521  std::string str = "[";
1522  for(size_t row = 0; row < val.matrix().rows() && row < maxRows; row++) {
1523  str += (row > 0 ? ";[" : "[");
1524  for(size_t col = 0; col < val.matrix().cols() && col < maxCols; col++) {
1525  if(col > 0) str += ",";
1526  str += numberToString<typename Derived::Scalar>(val.matrix()(row, col), precision);
1527  }
1528  str += (val.matrix().cols() > maxCols ? "...]" : "]");
1529  }
1530  str += (val.matrix().rows() > maxRows ? "...]" : "]");
1531  return str;
1532  }
1533  }
1534 # endif // #ifdef EIGENLAB_DEBUG
1535 #endif // #ifdef DEBUG
1536 
1537  template <typename Derived>
1538  std::string Parser<Derived>::trim(const std::string & str)
1539  {
1540  if(str.size() == 0) return str;
1541  std::string::const_iterator first = str.begin();
1542  std::string::const_iterator last = str.end() - 1;
1543  while(first < last && isspace(*first)) first++;
1544  while(first < last && isspace(*last)) last--;
1545  return std::string(first, last + 1);
1546  }
1547 
1548  template <typename Derived>
1549  std::vector<std::string> Parser<Derived>::split(const std::string & str, const char delimeter)
1550  {
1551  std::vector<std::string> args;
1552  std::string::const_iterator i0 = str.begin();
1553  for(std::string::const_iterator it = str.begin(); it != str.end(); it++) {
1554  if((*it) == delimeter) {
1555  args.push_back(trim(std::string(i0, it)));
1556  i0 = it + 1;
1557  }
1558  }
1559  args.push_back(std::string(i0, str.end()));
1560  return args;
1561  }
1562 
1563  template <typename Derived>
1564  template <typename T>
1565  bool Parser<Derived>::isNumber(const std::string & str, T * num)
1566  {
1567  std::istringstream iss(str);
1568  if(num)
1569  iss >> (* num);
1570  else {
1571  T number;
1572  iss >> number;
1573  }
1574  return (!iss.fail() && !iss.bad() && iss.eof());
1575  }
1576 
1577  template <typename Derived>
1578  template<typename T>
1579  T Parser<Derived>::stringToNumber(const std::string & str)
1580  {
1581  std::istringstream iss(str);
1582  T number;
1583  iss >> number;
1584  if(iss.fail() || iss.bad() || !iss.eof())
1585  throw std::runtime_error("Failed to convert " + str + " to a number.");
1586  return number;
1587  }
1588 
1589  template <typename Derived>
1590  template<typename T>
1591  std::string Parser<Derived>::numberToString(T num, int precision)
1592  {
1593  std::ostringstream oss;
1594  if(precision)
1595  oss << std::setprecision(precision) << num;
1596  else
1597  oss << num;
1598  return oss.str();
1599  }
1600 
1601 #ifdef DEBUG
1602  template <typename Derived>
1603  void
1604  Parser<Derived>::test_w_lt(size_t & numFails,
1605  typename Derived::Scalar & /* s */,
1606  Derived & a34,
1607  Derived & b34,
1608  Derived & /* c43 */,
1609  Derived & /* v */, std::true_type)
1610  {
1611  //
1612  // tests that only work if Derived::Scalar has operator<
1613  //
1614  Value<Derived> resultValue;
1615  Derived resultMatrix;
1616  Derived temp;
1617  std::cout << "Test min(a): ";
1618  resultValue = eval("min(a)");
1619  resultMatrix.setConstant(1, 1, a34.minCoeff());
1620  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1621  else { std::cout << "FAIL" << std::endl; ++numFails; }
1622 
1623  std::cout << "Test min(a, 0): ";
1624  resultValue = eval("min(a, 0)");
1625  resultMatrix = a34.colwise().minCoeff();
1626  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1627  else { std::cout << "FAIL" << std::endl; ++numFails; }
1628 
1629  std::cout << "Test min(a, 1): ";
1630  resultValue = eval("min(a, 1)");
1631  resultMatrix = a34.rowwise().minCoeff();
1632  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1633  else { std::cout << "FAIL" << std::endl; ++numFails; }
1634 
1635  std::cout << "Test max(a): ";
1636  resultValue = eval("max(a)");
1637  resultMatrix.setConstant(1, 1, a34.maxCoeff());
1638  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1639  else { std::cout << "FAIL" << std::endl; ++numFails; }
1640 
1641  std::cout << "Test max(a, 0): ";
1642  resultValue = eval("max(a, 0)");
1643  resultMatrix = a34.colwise().maxCoeff();
1644  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1645  else { std::cout << "FAIL" << std::endl; ++numFails; }
1646 
1647  std::cout << "Test max(a, 1): ";
1648  resultValue = eval("max(a, 1)");
1649  resultMatrix = a34.rowwise().maxCoeff();
1650  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1651  else { std::cout << "FAIL" << std::endl; ++numFails; }
1652 
1653  std::cout << "Test absmax(a): ";
1654  resultValue = eval("absmax(a)");
1655  resultMatrix.setConstant(1, 1, std::abs(a34.maxCoeff()) >= std::abs(a34.minCoeff()) ? a34.maxCoeff() : a34.minCoeff());
1656  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1657  else { std::cout << "FAIL" << std::endl; ++numFails; }
1658 
1659  std::cout << "Test absmax(a, 0): ";
1660  resultValue = eval("absmax(a, 0)");
1661  resultMatrix = a34.colwise().maxCoeff();
1662  temp = a34.colwise().minCoeff();
1663  for(Index i = 0; i < resultMatrix.size(); ++i) {
1664  if(std::abs(resultMatrix(i)) < std::abs(temp(i)))
1665  resultMatrix(i) = temp(i);
1666  }
1667  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1668  else { std::cout << "FAIL" << std::endl; ++numFails; }
1669 
1670  std::cout << "Test absmax(a, 1): ";
1671  resultValue = eval("absmax(a, 1)");
1672  resultMatrix = a34.rowwise().maxCoeff();
1673  temp = a34.rowwise().minCoeff();
1674  for(Index i = 0; i < resultMatrix.size(); ++i) {
1675  if(std::abs(resultMatrix(i)) < std::abs(temp(i)))
1676  resultMatrix(i) = temp(i);
1677  }
1678  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1679  else { std::cout << "FAIL" << std::endl; ++numFails; }
1680 
1681  std::cout << "Test cwiseMin(a, b): ";
1682  resultValue = eval("cwiseMin(a, b)");
1683  resultMatrix = a34.cwiseMin(b34);
1684  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1685  else { std::cout << "FAIL" << std::endl; ++numFails; }
1686 
1687  std::cout << "Test cwiseMax(a, b): ";
1688  resultValue = eval("cwiseMax(a, b)");
1689  resultMatrix = a34.cwiseMax(b34);
1690  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1691  else { std::cout << "FAIL" << std::endl; ++numFails; }
1692  }
1693 
1694  template <typename Derived>
1695  void
1696  Parser<Derived>::test_w_lt(size_t & /* numFails */,
1697  typename Derived::Scalar & /* s */,
1698  Derived & /* a34 */,
1699  Derived & /* b34 */,
1700  Derived & /* c43 */,
1701  Derived & /* v */, std::false_type)
1702  {
1703  // do nothing
1704  }
1705 
1706  template <typename Derived>
1707  size_t Parser<Derived>::test()
1708  {
1709  std::cout << std::endl;
1710  std::cout << "BEGIN unit test for EigenLab..." << std::endl;
1711  std::cout << "Make sure this function completes successfuly and prints the message "
1712  "'Successfully completed unit test for EigenLab with no failures.'" << std::endl;
1713  std::cout << std::endl;
1714 
1715  size_t numFails = 0;
1716  Value<Derived> resultValue;
1717  Derived resultMatrix;
1718  Derived temp;
1719  typename Derived::Scalar s = 2;
1720 
1721  Derived a34 = Derived::Random(3, 4);
1722  Derived b34 = Derived::Random(3, 4);
1723  Derived c43 = Derived::Random(4, 3);
1724  Derived v = Derived::Random(1, 10);
1725  //std::cout << "a34=" << std::endl << a34 << std::endl << std::endl;
1726  //std::cout << "b34=" << std::endl << b34 << std::endl << std::endl;
1727  //std::cout << "c43=" << std::endl << c43 << std::endl << std::endl;
1728  //std::cout << "v=" << std::endl << v << std::endl << std::endl;
1729 
1730  var("a").setShared(a34);
1731  var("b").setShared(b34);
1732  var("c").setShared(c43);
1733  var("v").setShared(v);
1734  var("s").setShared(& s);
1735 
1737  std::cout << "Testing basic operations..." << std::endl << std::endl;
1739 
1740  std::cout << "Test matrix addition a + b: ";
1741  resultValue = eval("a + b");
1742  resultMatrix = a34 + b34;
1743  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1744  else { std::cout << "FAIL" << std::endl; ++numFails; }
1745 
1746  std::cout << "Test matrix/scalar addition a + s: ";
1747  resultValue = eval("a + s");
1748  resultMatrix = a34.array() + s;
1749  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1750  else { std::cout << "FAIL" << std::endl; ++numFails; }
1751 
1752  std::cout << "Test scalar/matrix addition s + b: ";
1753  resultValue = eval("s + b");
1754  resultMatrix = b34.array() + s;
1755  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1756  else { std::cout << "FAIL" << std::endl; ++numFails; }
1757 
1758  std::cout << "Test matrix addition a .+ b: ";
1759  resultValue = eval("a .+ b");
1760  resultMatrix = a34 + b34;
1761  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1762  else { std::cout << "FAIL" << std::endl; ++numFails; }
1763 
1764  std::cout << "Test matrix/scalar addition a .+ s: ";
1765  resultValue = eval("a .+ s");
1766  resultMatrix = a34.array() + s;
1767  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1768  else { std::cout << "FAIL" << std::endl; ++numFails; }
1769 
1770  std::cout << "Test scalar/matrix addition s .+ b: ";
1771  resultValue = eval("s .+ b");
1772  resultMatrix = b34.array() + s;
1773  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1774  else { std::cout << "FAIL" << std::endl; ++numFails; }
1775 
1776  std::cout << "Test matrix subtraction a - b: ";
1777  resultValue = eval("a - b");
1778  resultMatrix = a34 - b34;
1779  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1780  else { std::cout << "FAIL" << std::endl; ++numFails; }
1781 
1782  std::cout << "Test matrix/scalar subtraction a - s: ";
1783  resultValue = eval("a - s");
1784  resultMatrix = a34.array() - s;
1785  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1786  else { std::cout << "FAIL" << std::endl; ++numFails; }
1787 
1788  std::cout << "Test scalar/matrix subtraction s - b: ";
1789  resultValue = eval("s - b");
1790  resultMatrix = (-b34.array()) + s;
1791  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1792  else { std::cout << "FAIL" << std::endl; ++numFails; }
1793 
1794  std::cout << "Test matrix subtraction a .- b: ";
1795  resultValue = eval("a .- b");
1796  resultMatrix = a34 - b34;
1797  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1798  else { std::cout << "FAIL" << std::endl; ++numFails; }
1799 
1800  std::cout << "Test matrix/scalar subtraction a .- s: ";
1801  resultValue = eval("a .- s");
1802  resultMatrix = a34.array() - s;
1803  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1804  else { std::cout << "FAIL" << std::endl; ++numFails; }
1805 
1806  std::cout << "Test scalar/matrix subtraction s .- b: ";
1807  resultValue = eval("s .- b");
1808  resultMatrix = (-b34.array()) + s;
1809  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1810  else { std::cout << "FAIL" << std::endl; ++numFails; }
1811 
1812  std::cout << "Test matrix negation -a: ";
1813  resultValue = eval("-a");
1814  resultMatrix = -a34;
1815  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1816  else { std::cout << "FAIL" << std::endl; ++numFails; }
1817 
1818  std::cout << "Test scalar negation -s: ";
1819  resultValue = eval("-s");
1820  resultMatrix.setConstant(1, 1, -s);
1821  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1822  else { std::cout << "FAIL" << std::endl; ++numFails; }
1823 
1824  std::cout << "Test matrix coefficient-wise multiplication a .* b: ";
1825  resultValue = eval("a .* b");
1826  resultMatrix = a34.array() * b34.array();
1827  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1828  else { std::cout << "FAIL" << std::endl; ++numFails; }
1829 
1830  std::cout << "Test matrix/scalar coefficient-wise multiplication a * s: ";
1831  resultValue = eval("a * s");
1832  resultMatrix = a34.array() * s;
1833  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1834  else { std::cout << "FAIL" << std::endl; ++numFails; }
1835 
1836  std::cout << "Test scalar/matrix coefficient-wise multiplication s * b: ";
1837  resultValue = eval("s * b");
1838  resultMatrix = b34.array() * s;
1839  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1840  else { std::cout << "FAIL" << std::endl; ++numFails; }
1841 
1842  std::cout << "Test matrix/scalar coefficient-wise multiplication a .* s: ";
1843  resultValue = eval("a .* s");
1844  resultMatrix = a34.array() * s;
1845  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1846  else { std::cout << "FAIL" << std::endl; ++numFails; }
1847 
1848  std::cout << "Test scalar/matrix coefficient-wise multiplication s .* b: ";
1849  resultValue = eval("s .* b");
1850  resultMatrix = b34.array() * s;
1851  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1852  else { std::cout << "FAIL" << std::endl; ++numFails; }
1853 
1854  std::cout << "Test matrix coefficient-wise division a ./ b: ";
1855  resultValue = eval("a ./ b");
1856  resultMatrix = a34.array() / b34.array();
1857  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1858  else { std::cout << "FAIL" << std::endl; ++numFails; }
1859 
1860  std::cout << "Test matrix/scalar coefficient-wise division a / s: ";
1861  resultValue = eval("a / s");
1862  resultMatrix = a34.array() / s;
1863  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1864  else { std::cout << "FAIL" << std::endl; ++numFails; }
1865 
1866  std::cout << "Test scalar/matrix coefficient-wise division s / b: ";
1867  resultValue = eval("s / b");
1868  resultMatrix = Derived::Constant(b34.rows(), b34.cols(), s).array() / b34.array();
1869  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1870  else { std::cout << "FAIL" << std::endl; ++numFails; }
1871 
1872  std::cout << "Test matrix/scalar coefficient-wise division a ./ s: ";
1873  resultValue = eval("a ./ s");
1874  resultMatrix = a34.array() / s;
1875  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1876  else { std::cout << "FAIL" << std::endl; ++numFails; }
1877 
1878  std::cout << "Test scalar/matrix coefficient-wise division s ./ b: ";
1879  resultValue = eval("s ./ b");
1880  resultMatrix = Derived::Constant(b34.rows(), b34.cols(), s).array() / b34.array();
1881  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1882  else { std::cout << "FAIL" << std::endl; ++numFails; }
1883 
1884  std::cout << "Test matrix coefficient-wise power a .^ b: ";
1885  resultValue = eval("abs(a) .^ b");
1886  resultMatrix = a34;
1887  for(Index i = 0; i < a34.size(); ++i)
1888  resultMatrix(i) = pow(std::abs(a34(i)), b34(i));
1889  // std::cout << std::endl;
1890  // std::cout << "a=" << std::endl << a34 << std::endl << std::endl;
1891  // std::cout << "b=" << std::endl << b34 << std::endl << std::endl;
1892  // std::cout << "val=" << std::endl << resultValue.matrix() << std::endl << std::endl;
1893  // std::cout << "mat=" << std::endl << resultMatrix << std::endl << std::endl;
1894  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1895  else { std::cout << "FAIL" << std::endl; ++numFails; }
1896 
1897  std::cout << "Test matrix/scalar coefficient-wise power a ^ s: ";
1898  resultValue = eval("abs(a) ^ s");
1899  resultMatrix = a34;
1900  for(Index i = 0; i < a34.size(); ++i)
1901  resultMatrix(i) = pow(std::abs(a34(i)), s);
1902  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1903  else { std::cout << "FAIL" << std::endl; ++numFails; }
1904 
1905  std::cout << "Test scalar/matrix coefficient-wise power s ^ b: ";
1906  resultValue = eval("s ^ b");
1907  resultMatrix = b34;
1908  for(Index i = 0; i < b34.size(); ++i)
1909  resultMatrix(i) = pow(s, b34(i));
1910  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1911  else { std::cout << "FAIL" << std::endl; ++numFails; }
1912 
1913  std::cout << "Test matrix/scalar coefficient-wise power a .^ s: ";
1914  resultValue = eval("abs(a) .^ s");
1915  resultMatrix = a34;
1916  for(Index i = 0; i < a34.size(); ++i)
1917  resultMatrix(i) = pow(std::abs(a34(i)), s);
1918  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1919  else { std::cout << "FAIL" << std::endl; ++numFails; }
1920 
1921  std::cout << "Test scalar/matrix coefficient-wise power s .^ b: ";
1922  resultValue = eval("s .^ b");
1923  resultMatrix = b34;
1924  for(Index i = 0; i < b34.size(); ++i)
1925  resultMatrix(i) = pow(s, b34(i));
1926  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1927  else { std::cout << "FAIL" << std::endl; ++numFails; }
1928 
1929  std::cout << "Test matrix multiplication a * b: ";
1930  resultValue = eval("a * c");
1931  resultMatrix = a34 * c43;
1932  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1933  else { std::cout << "FAIL" << std::endl; ++numFails; }
1934 
1935  std::cout << "Test grouped subexpression (a + b) * c: ";
1936  resultValue = eval("(a + b) * c");
1937  resultMatrix = (a34 + b34) * c43;
1938  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1939  else { std::cout << "FAIL" << std::endl; ++numFails; }
1940 
1941  std::cout << "Test nested groups ((a + (b * 3 + 1)) * c).^2: ";
1942  resultValue = eval("((a + (b * 3 + 1)) * c).^2");
1943  resultMatrix = ((a34.array() + (b34.array() * 3 + 1)).matrix() * c43).array().pow(2);
1944  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1945  else { std::cout << "FAIL" << std::endl; ++numFails; }
1946 
1948  std::cout << std::endl << "Testing coefficient and submatrix block access..." << std::endl << std::endl;
1950 
1951  std::cout << "Test matrix coefficient access a(i,j): ";
1952  resultValue = eval("a(1,2)");
1953  resultMatrix.setConstant(1, 1, a34(1,2));
1954  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1955  else { std::cout << "FAIL" << std::endl; ++numFails; }
1956 
1957  std::cout << "Test submatrix block access a(i:p,j:q): ";
1958  resultValue = eval("a(1:2,2:3)");
1959  resultMatrix = a34.block(1, 2, 2, 2);
1960  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1961  else { std::cout << "FAIL" << std::endl; ++numFails; }
1962 
1963  std::cout << "Test submatrix block access using 'end' and ':' identifiers a(i:end,:): ";
1964  resultValue = eval("a(1:end,:)");
1965  resultMatrix = a34.block(1, 0, a34.rows() - 1, a34.cols());
1966  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1967  else { std::cout << "FAIL" << std::endl; ++numFails; }
1968 
1969  std::cout << "Test submatrix block access using subexpressions: ";
1970  resultValue = eval("a(2-1:2-1,0+1:3-1)");
1971  resultMatrix = a34.block(1, 1, 1, 2);
1972  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1973  else { std::cout << "FAIL" << std::endl; ++numFails; }
1974 
1975  std::cout << "Test submatrix block access using subexpressions with 'end' keyword: ";
1976  resultValue = eval("a(2-1:end-1,0+1:end-1)");
1977  resultMatrix = a34.block(1, 1, a34.rows() - 2, a34.cols() - 2);
1978  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1979  else { std::cout << "FAIL" << std::endl; ++numFails; }
1980 
1981  std::cout << "Test vector coefficient access v(i): ";
1982  resultValue = eval("v(5)");
1983  resultMatrix.setConstant(1, 1, v(5));
1984  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1985  else { std::cout << "FAIL" << std::endl; ++numFails; }
1986 
1987  std::cout << "Test subvector segment access v(i:j): ";
1988  resultValue = eval("v(3:6)");
1989  resultMatrix = v.block(0, 3, 1, 4);
1990  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1991  else { std::cout << "FAIL" << std::endl; ++numFails; }
1992 
1993  std::cout << "Test subvector segment access using 'end' identifier v(i:end): ";
1994  resultValue = eval("v(5:end)");
1995  resultMatrix = v.block(0, 5, 1, v.cols() - 5);
1996  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
1997  else { std::cout << "FAIL" << std::endl; ++numFails; }
1998 
1999  std::cout << "Test subvector segment access using ':' identifier v(:): ";
2000  resultValue = eval("v(:)");
2001  resultMatrix = v;
2002  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2003  else { std::cout << "FAIL" << std::endl; ++numFails; }
2004 
2005  std::cout << "Test subvector segment access using subexpressions: ";
2006  resultValue = eval("v(3-1:5+2)");
2007  resultMatrix = v.block(0, 2, 1, 6);
2008  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2009  else { std::cout << "FAIL" << std::endl; ++numFails; }
2010 
2011  std::cout << "Test subvector segment access using subexpressions with 'end' keyword: ";
2012  resultValue = eval("v((end-8)*2:end-3)");
2013  resultMatrix = v.block(0, 2, 1, 5);
2014  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2015  else { std::cout << "FAIL" << std::endl; ++numFails; }
2016 
2018  std::cout << std::endl << "Testing vector/matrix expressions..." << std::endl << std::endl;
2020 
2021  std::cout << "Test numeric range expression [i:j]: ";
2022  resultValue = eval("[2:5]");
2023  resultMatrix.resize(1, 4);
2024  resultMatrix << 2, 3, 4, 5;
2025  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2026  else { std::cout << "FAIL" << std::endl; ++numFails; }
2027 
2028  std::cout << "Test numeric range expression [i:s:j]: ";
2029  resultValue = eval("[2:2:10]");
2030  resultMatrix.resize(1, 5);
2031  resultMatrix << 2, 4, 6, 8, 10;
2032  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2033  else { std::cout << "FAIL" << std::endl; ++numFails; }
2034 
2035  std::cout << "Test numeric range with subexpressions: ";
2036  resultValue = eval("[6-2:5*2-3]");
2037  std::cout << "val=" << std::endl << resultValue.matrix() << std::endl << std::endl;
2038  resultMatrix.resize(1, 4);
2039  resultMatrix << 4, 5, 6, 7;
2040  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2041  else { std::cout << "FAIL" << std::endl; ++numFails; }
2042 
2043  std::cout << "Test matrix expression [[1, 2]; [3, 4]]: ";
2044  resultValue = eval("[[1, 2]; [3, 4]]");
2045  resultMatrix.resize(2, 2);
2046  resultMatrix << 1, 2, 3, 4;
2047  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2048  else { std::cout << "FAIL" << std::endl; ++numFails; }
2049 
2050  std::cout << "Test matrix expression [ 1, 2, 3; 4:6 ]: ";
2051  resultValue = eval("[1, 2, 3; 4:6]");
2052  resultMatrix.resize(2, 3);
2053  resultMatrix << 1, 2, 3, 4, 5, 6;
2054  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2055  else { std::cout << "FAIL" << std::endl; ++numFails; }
2056 
2058  std::cout << std::endl << "Testing coefficient-wise functions..." << std::endl << std::endl;
2060 
2061  std::cout << "Test coefficient-wise abs(a): ";
2062  resultValue = eval("abs(a)");
2063  resultMatrix.resize(3, 4);
2064  resultMatrix.real() = a34.array().abs();
2065  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2066  else { std::cout << "FAIL" << std::endl; ++numFails; }
2067 
2068  std::cout << "Test coefficient-wise sqrt(a): ";
2069  resultValue = eval("sqrt(abs(a))");
2070  resultMatrix.real() = a34.array().abs().sqrt();
2071  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2072  else { std::cout << "FAIL" << std::endl; ++numFails; }
2073 
2074  std::cout << "Test coefficient-wise exp(a): ";
2075  resultValue = eval("exp(abs(a) + 0.001)");
2076  resultMatrix.real() = (a34.array().abs() + 0.001).exp();
2077  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2078  else { std::cout << "FAIL" << std::endl; ++numFails; }
2079 
2080  std::cout << "Test coefficient-wise log(a): ";
2081  resultValue = eval("log(abs(a) + 0.001)");
2082  resultMatrix.real() = (a34.array().abs() + 0.001).log();
2083  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2084  else { std::cout << "FAIL" << std::endl; ++numFails; }
2085 
2086  std::cout << "Test coefficient-wise log10(a): ";
2087  resultValue = eval("log10(abs(a) + 0.001)");
2088  resultMatrix.real() = (a34.array().abs() + 0.001).log() * (1.0 / log(10));
2089  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2090  else { std::cout << "FAIL" << std::endl; ++numFails; }
2091 
2092  std::cout << "Test coefficient-wise sin(a): ";
2093  resultValue = eval("sin(a)");
2094  resultMatrix = a34.array().sin();
2095  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2096  else { std::cout << "FAIL" << std::endl; ++numFails; }
2097 
2098  std::cout << "Test coefficient-wise cos(a): ";
2099  resultValue = eval("cos(a)");
2100  resultMatrix = a34.array().cos();
2101  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2102  else { std::cout << "FAIL" << std::endl; ++numFails; }
2103 
2104  std::cout << "Test coefficient-wise tan(a): ";
2105  resultValue = eval("tan(a)");
2106  resultMatrix = a34.array().tan();
2107  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2108  else { std::cout << "FAIL" << std::endl; ++numFails; }
2109 
2110  std::cout << "Test coefficient-wise asin(a): ";
2111  resultValue = eval("asin(a)");
2112  resultMatrix = a34.array().asin();
2113  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2114  else { std::cout << "FAIL" << std::endl; ++numFails; }
2115 
2116  std::cout << "Test coefficient-wise acos(a): ";
2117  resultValue = eval("acos(a)");
2118  resultMatrix = a34.array().acos();
2119  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2120  else { std::cout << "FAIL" << std::endl; ++numFails; }
2121 
2123  std::cout << std::endl << "Testing matrix reduction functions..." << std::endl << std::endl;
2125 
2126  std::cout << "Test trace(a): ";
2127  resultValue = eval("trace(a)");
2128  resultMatrix.setConstant(1, 1, a34.trace());
2129  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2130  else { std::cout << "FAIL" << std::endl; ++numFails; }
2131 
2132  std::cout << "Test norm(a): ";
2133  resultValue = eval("norm(a)");
2134  resultMatrix.setConstant(1, 1, a34.norm());
2135  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2136  else { std::cout << "FAIL" << std::endl; ++numFails; }
2137 
2138  std::cout << "Test size(a, 0): ";
2139  resultValue = eval("size(a, 0)");
2140  resultMatrix.setConstant(1, 1, a34.rows());
2141  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2142  else { std::cout << "FAIL" << std::endl; ++numFails; }
2143 
2144  std::cout << "Test size(a, 1): ";
2145  resultValue = eval("size(a, 1)");
2146  resultMatrix.setConstant(1, 1, a34.cols());
2147  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2148  else { std::cout << "FAIL" << std::endl; ++numFails; }
2149 
2150  test_w_lt(numFails, s, a34, b34, c43, v, has_operator_lt<typename Derived::Scalar>());
2151 
2152  std::cout << "Test mean(a): ";
2153  resultValue = eval("mean(a)");
2154  resultMatrix.setConstant(1, 1, a34.mean());
2155  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2156  else { std::cout << "FAIL" << std::endl; ++numFails; }
2157 
2158  std::cout << "Test mean(a, 0): ";
2159  resultValue = eval("mean(a, 0)");
2160  resultMatrix = a34.colwise().mean();
2161  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2162  else { std::cout << "FAIL" << std::endl; ++numFails; }
2163 
2164  std::cout << "Test mean(a, 1): ";
2165  resultValue = eval("mean(a, 1)");
2166  resultMatrix = a34.rowwise().mean();
2167  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2168  else { std::cout << "FAIL" << std::endl; ++numFails; }
2169 
2170  std::cout << "Test sum(a): ";
2171  resultValue = eval("sum(a)");
2172  resultMatrix.setConstant(1, 1, a34.sum());
2173  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2174  else { std::cout << "FAIL" << std::endl; ++numFails; }
2175 
2176  std::cout << "Test sum(a, 0): ";
2177  resultValue = eval("sum(a, 0)");
2178  resultMatrix = a34.colwise().sum();
2179  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2180  else { std::cout << "FAIL" << std::endl; ++numFails; }
2181 
2182  std::cout << "Test sum(a, 1): ";
2183  resultValue = eval("sum(a, 1)");
2184  resultMatrix = a34.rowwise().sum();
2185  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2186  else { std::cout << "FAIL" << std::endl; ++numFails; }
2187 
2188  std::cout << "Test prod(a): ";
2189  resultValue = eval("prod(a)");
2190  resultMatrix.setConstant(1, 1, a34.prod());
2191  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2192  else { std::cout << "FAIL" << std::endl; ++numFails; }
2193 
2194  std::cout << "Test prod(a, 0): ";
2195  resultValue = eval("prod(a, 0)");
2196  resultMatrix = a34.colwise().prod();
2197  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2198  else { std::cout << "FAIL" << std::endl; ++numFails; }
2199 
2200  std::cout << "Test prod(a, 1): ";
2201  resultValue = eval("prod(a, 1)");
2202  resultMatrix = a34.rowwise().prod();
2203  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2204  else { std::cout << "FAIL" << std::endl; ++numFails; }
2205 
2207  std::cout << std::endl << "Testing matrix functions..." << std::endl << std::endl;
2209 
2210  std::cout << "Test transpose(a): ";
2211  resultValue = eval("transpose(a)");
2212  resultMatrix = a34.transpose();
2213  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2214  else { std::cout << "FAIL" << std::endl; ++numFails; }
2215 
2216  std::cout << "Test conjugate(a): ";
2217  resultValue = eval("conjugate(a)");
2218  resultMatrix = a34.conjugate();
2219  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2220  else { std::cout << "FAIL" << std::endl; ++numFails; }
2221 
2222  std::cout << "Test adjoint(a): ";
2223  resultValue = eval("adjoint(a)");
2224  resultMatrix = a34.adjoint();
2225  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2226  else { std::cout << "FAIL" << std::endl; ++numFails; }
2227 
2229  std::cout << std::endl << "Testing matrix initializers..." << std::endl << std::endl;
2231 
2232  std::cout << "Test zeros(5): ";
2233  resultValue = eval("zeros(5)");
2234  resultMatrix = Derived::Zero(5, 5);
2235  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2236  else { std::cout << "FAIL" << std::endl; ++numFails; }
2237 
2238  std::cout << "Test ones(5): ";
2239  resultValue = eval("ones(5)");
2240  resultMatrix = Derived::Ones(5, 5);
2241  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2242  else { std::cout << "FAIL" << std::endl; ++numFails; }
2243 
2244  std::cout << "Test eye(5): ";
2245  resultValue = eval("eye(5)");
2246  resultMatrix = Derived::Identity(5, 5);
2247  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2248  else { std::cout << "FAIL" << std::endl; ++numFails; }
2249 
2250  try {
2251  std::cout << "Test zeros(5.2): ";
2252  resultValue = eval("zeros(5.2)"); // <-- Should NOT succeed!!!
2253  std::cout << "FAIL" << std::endl; ++numFails;
2254  } catch(std::runtime_error &err) {
2255  std::cout << err.what() << std::endl;
2256  std::cout << "Exception caught, so we're OK" << std::endl;
2257  }
2258 
2259  try {
2260  std::cout << "Test eye(-3): ";
2261  resultValue = eval("eye(-3)"); // <-- Should NOT succeed!!!
2262  std::cout << "FAIL" << std::endl; ++numFails;
2263  } catch(std::runtime_error &err) {
2264  std::cout << err.what() << std::endl;
2265  std::cout << "Exception caught, so we're OK" << std::endl;
2266  }
2267 
2268  std::cout << "Test zeros(4,7): ";
2269  resultValue = eval("zeros(4,7)");
2270  resultMatrix = Derived::Zero(4, 7);
2271  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2272  else { std::cout << "FAIL" << std::endl; ++numFails; }
2273 
2274  std::cout << "Test ones(4,7): ";
2275  resultValue = eval("ones(4,7)");
2276  resultMatrix = Derived::Ones(4, 7);
2277  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2278  else { std::cout << "FAIL" << std::endl; ++numFails; }
2279 
2280  std::cout << "Test eye(4,7): ";
2281  resultValue = eval("eye(4,7)");
2282  resultMatrix = Derived::Identity(4, 7);
2283  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2284  else { std::cout << "FAIL" << std::endl; ++numFails; }
2285 
2287  std::cout << std::endl << "Testing variable assignment..." << std::endl << std::endl;
2289 
2290  std::cout << "Test assigning to a variable with the same dimensions a = b: ";
2291  eval("a = b");
2292  if(a34.isApprox(b34)) std::cout << "OK" << std::endl;
2293  else { std::cout << "FAIL" << std::endl; ++numFails; }
2294 
2295  std::cout << "Test assigning to a variable with different dimensions a = c: ";
2296  eval("a = c");
2297  if(var("a").matrix().isApprox(c43) && a34.isApprox(b34)) std::cout << "OK" << std::endl;
2298  else { std::cout << "FAIL" << std::endl; ++numFails; }
2299  var("a").setShared(a34);
2300 
2301  std::cout << "Test creating a new variable x = [1,2;3,4]: ";
2302  resultValue = eval("x = [1,2;3,4]");
2303  if(var("x").matrix().isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2304  else { std::cout << "FAIL" << std::endl; ++numFails; }
2305 
2306  std::cout << "Test assigning to a variable coefficient a(i,j) = s: ";
2307  eval("a(1, 2) = s");
2308  if(a34(1, 2) == s) std::cout << "OK" << std::endl;
2309  else { std::cout << "FAIL" << std::endl; ++numFails; }
2310 
2311  std::cout << "Test assigning to a variable submatrix block a(0:1,1:2) = x: ";
2312  resultValue = eval("a(0:1,1:2) = x");
2313  if(a34.block(0,1,2,2).isApprox(var("x").matrix())) std::cout << "OK" << std::endl;
2314  else { std::cout << "FAIL" << std::endl; ++numFails; }
2315 
2316  try {
2317  std::cout << "Test bad function call: ";
2318  resultValue = eval("foobar(-3)"); // <-- Should NOT succeed!!!
2319  std::cout << "FAIL" << std::endl; ++numFails;
2320  } catch(std::runtime_error &err) {
2321  std::cout << err.what() << std::endl;
2322  std::cout << "Exception caught, so we're OK" << std::endl;
2323  }
2325  std::cout << std::endl << "Testing fp parsing..." << std::endl << std::endl;
2327 
2328  std::cout << "Test assigning 1.2e-3: ";
2329  resultValue = eval("s = 1.2e-3");
2330  resultMatrix.setConstant(1, 1, typename Derived::Scalar(1.2e-3));
2331  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2332  else { std::cout << "FAIL" << std::endl; ++numFails; }
2333 
2334  std::cout << "Test assigning 1.e3: ";
2335  resultValue = eval("s = 1.e3");
2336  resultMatrix.setConstant(1, 1, typename Derived::Scalar(1000));
2337  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2338  else { std::cout << "FAIL" << std::endl; ++numFails; }
2339 
2340  std::cout << "Test assigning 12.34e05: ";
2341  resultValue = eval("s = 12.34e05");
2342  resultMatrix.setConstant(1, 1, typename Derived::Scalar(123.4e4));
2343  if(resultMatrix.isApprox(resultValue.matrix())) std::cout << "OK" << std::endl;
2344  else { std::cout << "FAIL" << std::endl; ++numFails; }
2345 
2346  std::cout << std::endl;
2347  if(numFails == 0)
2348  std::cout << "Successfully completed unit test for EigenLab with no failures." << std::endl;
2349  else
2350  std::cout << "Completed unit test for EigenLab with " << numFails << " failures (see above)." << std::endl;
2351  std::cout << std::endl;
2352  return numFails;
2353  }
2354 #endif // #ifdef DEBUG
2355 
2356 } // namespace EigenLab
2357 
2358 #endif // #ifndef EigenLab_H
void evalMatrixExpression(const std::string &str, Value< Derived > &mat)
Definition: EigenLab.h:581
std::map< std::string, Value< Derived > > ValueMap
Definition: EigenLab.h:135
void evalFunction(const std::string &name, std::vector< std::string > &args, Value< Derived > &result)
Definition: EigenLab.h:805
ValueMap mVariables
Definition: EigenLab.h:139
void setShared(const Value &val)
Definition: EigenLab.h:102
void setShared(const Derived &mat)
Definition: EigenLab.h:101
const Eigen::Map< Derived > & matrix() const
Definition: EigenLab.h:69
std::vector< std::string > mFunctions
Definition: EigenLab.h:143
static std::vector< std::string > split(const std::string &str, const char delimeter)
Definition: EigenLab.h:1549
void setShared(const typename Derived::Scalar *data, size_t rows=1, size_t cols=1)
Definition: EigenLab.h:100
std::string mOperators2
Definition: EigenLab.h:142
typename std::is_same< bool, decltype(test< T >(0))>::type type
Definition: EigenLab.h:120
Value< Eigen::MatrixXd > ValueXd
Definition: EigenLab.h:108
void clearCachedExpressions()
Definition: EigenLab.h:174
void setLocal(const Eigen::MatrixBase< Derived > &mat)
Definition: EigenLab.h:92
XmlRpcServer s
void setLocal(const typename Derived::Scalar *data, size_t rows=1, size_t cols=1)
Definition: EigenLab.h:94
Value< Derived > & var(const std::string &name)
Definition: EigenLab.h:163
Chunk(const std::string &str="", int t=-1, const Value< Derived > &val=Value< Derived >())
Definition: EigenLab.h:151
Parser< Eigen::MatrixXd > ParserXd
Definition: EigenLab.h:225
Derived mLocal
Definition: EigenLab.h:54
bool hasVar(const std::string &name)
Definition: EigenLab.h:166
Derived::Index Index
Definition: EigenLab.h:155
Value< Derived > value
Definition: EigenLab.h:149
void evalNumericRange(const std::string &str, Value< Derived > &mat)
Definition: EigenLab.h:1030
void mapLocal()
Definition: EigenLab.h:82
void evalIndexRange(const std::string &str, int *first, int *last, int numIndices)
Definition: EigenLab.h:525
Value(const Value &val)
Definition: EigenLab.h:105
bool evalFunction_2_lt(const std::string &name, Value< Derived > &arg0, Value< Derived > &arg1, Value< Derived > &result, std::false_type)
Definition: EigenLab.h:799
Eigen::Map< Derived > & matrix()
Definition: EigenLab.h:68
void evalPowers(ChunkArray &chunks)
Definition: EigenLab.h:1181
std::vector< std::string > splitArguments(const std::string &str, const char delimeter) const
Definition: EigenLab.h:508
static std::string numberToString(T num, int precision=0)
Definition: EigenLab.h:1591
Value< Eigen::MatrixXf > ValueXf
Definition: EigenLab.h:109
void clearVar(const std::string &name)
Definition: EigenLab.h:169
void setCacheExpressions(bool b)
Definition: EigenLab.h:173
Value(const typename Derived::Scalar s)
Definition: EigenLab.h:89
Derived & local()
Definition: EigenLab.h:75
bool evalFunction_1_lt(const std::string &name, Value< Derived > &arg, Value< Derived > &result, std::false_type)
Definition: EigenLab.h:703
void splitEquationIntoChunks(const std::string &expression, ChunkArray &chunks, std::string &code)
Definition: EigenLab.h:316
Parser< Eigen::MatrixXf > ParserXf
Definition: EigenLab.h:226
Value< Eigen::MatrixXi > ValueXi
Definition: EigenLab.h:110
Value & operator=(const typename Derived::Scalar s)
Definition: EigenLab.h:95
std::map< std::string, ChunkArray > mCachedChunkedExpressions
Definition: EigenLab.h:157
void evalAddition(ChunkArray &chunks)
Definition: EigenLab.h:1339
const Derived & local() const
Definition: EigenLab.h:76
unsigned int step
static T stringToNumber(const std::string &str)
Definition: EigenLab.h:1579
void copySharedToLocal()
Definition: EigenLab.h:85
Parser< Eigen::MatrixXi > ParserXi
Definition: EigenLab.h:227
static std::string trim(const std::string &str)
Definition: EigenLab.h:1538
bool mCacheChunkedExpressions
Definition: EigenLab.h:156
ValueMap & vars()
Definition: EigenLab.h:162
bool isOperator(const char c) const
Definition: EigenLab.h:196
void setLocal(const Value &val)
Definition: EigenLab.h:93
void evalIndices(ChunkArray &chunks)
Definition: EigenLab.h:1100
bool cacheExpressions() const
Definition: EigenLab.h:172
void setLocal(const typename Derived::Scalar s)
Definition: EigenLab.h:91
Value(const Derived &mat)
Definition: EigenLab.h:90
bool isVariable(const std::string &name) const
Definition: EigenLab.h:195
void evalAssignment(ChunkArray &chunks)
Definition: EigenLab.h:1421
void evalNegations(ChunkArray &chunks)
Definition: EigenLab.h:1139
std::string::const_iterator findClosingBracket(const std::string &str, const std::string::const_iterator openingBracket, const char closingBracket) const
Definition: EigenLab.h:495
void evalMultiplication(ChunkArray &chunks)
Definition: EigenLab.h:1248
bool isLocal() const
Definition: EigenLab.h:79
Value(const typename Derived::Scalar *data, size_t rows=1, size_t cols=1)
Definition: EigenLab.h:99
bool isFunction(const std::string &str) const
Definition: EigenLab.h:198
Eigen::Map< Derived > mShared
Definition: EigenLab.h:59
std::vector< Chunk > ChunkArray
Definition: EigenLab.h:154
Value< Derived > eval(const std::string &expression)
Definition: EigenLab.h:283
static bool isNumber(const std::string &str, T *num=0)
Definition: EigenLab.h:1565


grid_map_filters
Author(s): Péter Fankhauser , Martin Wermelinger
autogenerated on Tue Jun 25 2019 20:02:22