13 #include <Eigen/Dense> 15 #include <type_traits> 27 #ifndef EIGENLAB_DEBUG 49 template <
typename Derived = Eigen::MatrixXd>
82 inline void mapLocal() {
new (&
mShared) Eigen::Map<Derived>(mLocal.data(), mLocal.rows(), mLocal.cols()); mIsLocal =
true; }
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(); }
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; }
105 Value(
const Value & val) : mLocal(1, 1), mShared(mLocal.data(), mLocal.rows(), mLocal.cols()) { (* this) = val; }
117 static auto test(U*) -> decltype(std::declval<U>() < std::declval<U>());
119 static auto test(...) -> std::false_type;
120 using type =
typename std::is_same<bool, decltype(test<T>(0))>
::type;
130 template <
typename Derived = Eigen::MatrixXd>
135 typedef std::map<std::string, Value<Derived> >
ValueMap;
155 typedef typename Derived::Index
Index;
162 inline ValueMap &
vars() {
return mVariables; }
166 inline bool hasVar(
const std::string & name) {
return isVariable(name); }
169 inline void clearVar(
const std::string & name) {
typename ValueMap::iterator it = mVariables.find(name);
if(it != mVariables.end()) mVariables.erase(it); }
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);
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);
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);
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);
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);
232 template <
typename Derived>
234 mOperators1(
"+-*/^()[]="),
235 mOperators2(
".+.-.*./.^"),
236 mCacheChunkedExpressions(false)
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");
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");
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");
272 mFunctions.push_back(
"transpose");
273 mFunctions.push_back(
"conjugate");
274 mFunctions.push_back(
"adjoint");
277 mFunctions.push_back(
"zeros");
278 mFunctions.push_back(
"ones");
279 mFunctions.push_back(
"eye");
282 template <
typename Derived>
286 # ifdef EIGENLAB_DEBUG 287 std::cout <<
"---" << std::endl;
288 std::cout <<
"EXPRESSION: " << expression << std::endl;
300 if(chunks.size() != 1)
301 throw std::runtime_error(
"Failed to reduce expression '" + expression +
"' to a single value.");
303 # ifdef EIGENLAB_DEBUG 304 std::cout <<
"---" << std::endl;
309 throw std::runtime_error(
"Unknown variable '" + chunks[0].field +
"'.");
310 return mVariables[chunks[0].field];
312 return chunks[0].value;
315 template <
typename Derived>
322 # ifdef EIGENLAB_DEBUG 323 std::cout <<
"CACHED CHUNKS: "; printChunks(chunks); std::cout << std::endl;
330 for(std::string::const_iterator it = expression.begin(); it != expression.end();)
332 int prevType = (chunks.size() ? chunks.back().type : -1);
337 }
else if(ci ==
'(' && (prevType ==
VALUE || prevType ==
VARIABLE)) {
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);
345 throw std::runtime_error(
"Unknown variable '" + chunks.back().field +
"'.");
346 chunks.back().value.setShared(
var(chunks.back().field));
349 int rows = int(chunks.back().value.matrix().rows());
350 int cols = int(chunks.back().value.matrix().cols());
352 if(args.size() == 1) {
355 chunks.back().row0 = first;
356 chunks.back().col0 = 0;
357 chunks.back().rows = last + 1 - first;
358 chunks.back().cols = 1;
359 }
else if(rows == 1) {
361 chunks.back().row0 = 0;
362 chunks.back().col0 = first;
363 chunks.back().rows = 1;
364 chunks.back().cols = last + 1 - first;
366 throw std::runtime_error(
"Missing row or column indices for '(" + chunks.back().field +
"(" + field +
")'.");
368 }
else if(args.size() == 2) {
370 chunks.back().row0 = first;
371 chunks.back().rows = last + 1 - first;
373 chunks.back().col0 = first;
374 chunks.back().cols = last + 1 - first;
376 throw std::runtime_error(
"Invalid index expression '" + chunks.back().field +
"(" + field +
")'.");
380 }
else if(ci ==
'(' && prevType ==
FUNCTION) {
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);
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");
392 }
else if(ci ==
'(') {
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);
399 code += (chunks.back().value.matrix().size() == 1 ?
"n" :
"M");
401 }
else if(ci ==
'[') {
404 throw std::runtime_error(
"Invalid operation '" + chunks.back().field + std::string(1, ci) +
"'.");
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);
409 chunks.push_back(
Chunk(
"[" + field +
"]", prevType =
VALUE));
411 code += (chunks.back().value.matrix().size() == 1 ?
"n" :
"M");
413 }
else if(it + 1 != expression.end() &&
isOperator(std::string(it, it + 2))) {
415 std::string field = std::string(it, it + 2);
421 std::string field = std::string(1, ci);
427 std::string::const_iterator jt = it + 1;
429 unsigned char state = 1;
430 for(std::string::const_iterator kt = it; state && kt != expression.end(); kt++) {
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;
438 const static char nextstate[9][5] = {{0},
449 state = nextstate[state][token];
450 if (state == 8) jt = kt;
452 for(; jt != expression.end(); jt++) {
456 std::string field =
trim(std::string(it, jt));
458 throw std::runtime_error(
"Invalid operation '" + chunks.back().field + field +
"'.");
460 if(field.find(
":") != std::string::npos) {
462 chunks.push_back(
Chunk(field, prevType =
VALUE));
464 code += (chunks.back().value.matrix().size() == 1 ?
"n" :
"M");
465 }
else if(isNumber<double>(field, & num)) {
472 code += (mVariables[field].matrix().size() == 1 ?
"vn" :
"vM");
485 # ifdef EIGENLAB_DEBUG 486 std::cout <<
"CHUNKS: "; printChunks(chunks); std::cout << std::endl;
487 std::cout <<
"CODE: " << code << std::endl;
494 template <
typename Derived>
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;
507 template <
typename Derived>
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++) {
515 else if((*it) == delimeter) {
516 args.push_back(
trim(std::string(i0, it)));
520 args.push_back(std::string(i0, str.end()));
524 template <
typename Derived>
528 throw std::runtime_error(
"Empty index range.");
532 for(std::string::const_iterator it = str.begin(); it != str.end(); 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()) {
538 (* last) = numIndices - 1;
541 if(firstStr.empty() || lastStr.empty())
542 throw std::runtime_error(
"Missing indices for '" + str +
"'.");
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);
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);
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);
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);
566 std::string firstStr = str;
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);
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);
580 template <
typename Derived>
585 std::vector<std::vector<typename Derived::Scalar> > temp;
587 size_t row0 = 0, col0 = 0, nrows = 0, ncols = 0;
589 for(
size_t i = 0; i < rows.size(); i++) {
591 if(rows[i][0] ==
'[' && rows[i].back() ==
']') rows[i] = rows[i].substr(1,
int(rows[i].size()) - 2);
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] +
"'.");
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);
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;
625 throw std::runtime_error(
"Invalid matrix definition '[" + str +
"]'. Invalid range '" + cols[j] +
"'.");
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);
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;
645 throw std::runtime_error(
"Invalid matrix definition '[" + str +
"]'. Invalid range '" + cols[j] +
"'.");
649 submatrix =
eval(cols[j]);
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);
661 col0 += submatrix.
matrix().cols();
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.");
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];
678 template <
typename Derived>
684 }
else if(name ==
"minOfFinites") {
687 }
else if(name ==
"max") {
690 }
else if(name ==
"maxOfFinites") {
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);
702 template <
typename Derived>
708 template <
typename Derived>
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 +
"(...)'.");
718 result.
local() = arg0.
matrix().colwise().minCoeff();
721 }
else if(dim == 1) {
722 result.
local() = arg0.
matrix().rowwise().minCoeff();
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 +
"(...)'.");
733 result.
local() = arg0.
matrix().colwise().maxCoeff();
736 }
else if(dim == 1) {
737 result.
local() = arg0.
matrix().rowwise().maxCoeff();
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 +
"(...)'.");
748 result.
local() = arg0.
matrix().colwise().maxCoeff();
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);
756 }
else if(dim == 1) {
757 result.
local() = arg0.
matrix().rowwise().maxCoeff();
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);
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);
778 throw std::runtime_error(
"Invalid dimension argument for function '" + name +
"(...)'.");
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);
792 throw std::runtime_error(
"Invalid dimension argument for function '" + name +
"(...)'.");
798 template <
typename Derived>
804 template <
typename Derived>
807 if(args.size() == 1) {
810 result.
local() = arg.
matrix().array().abs().template cast<typename Derived::Scalar>();
813 }
else if(name ==
"sqrt") {
817 }
else if(name ==
"square") {
821 }
else if(name ==
"exp") {
825 }
else if(name ==
"log") {
829 }
else if(name ==
"log10") {
831 result.
local() *= (1.0 / log(10));
834 }
else if(name ==
"sin") {
838 }
else if(name ==
"cos") {
842 }
else if(name ==
"tan") {
846 }
else if(name ==
"acos") {
850 }
else if(name ==
"asin") {
854 }
else if(name ==
"trace") {
857 }
else if(name ==
"norm") {
862 }
else if(name ==
"mean") {
865 }
else if(name ==
"meanOfFinites") {
868 }
else if(name ==
"sum") {
871 }
else if(name ==
"sumOfFinites") {
874 }
else if(name ==
"prod") {
877 }
else if(name ==
"numberOfFinites") {
880 }
else if(name ==
"transpose") {
884 }
else if(name ==
"conjugate") {
888 }
else if(name ==
"adjoint") {
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);
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);
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);
921 throw std::runtime_error(
"Invalid function '" + name +
"(" + args[0] +
")'.");
923 }
else if(args.size() == 2) {
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] +
")'.");
935 }
else if(dim == 1) {
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] +
")'.");
951 }
else if(dim == 1) {
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] +
")'.");
966 }
else if(dim == 1) {
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] +
")'.");
981 }
else if(dim == 1) {
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);
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);
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);
1017 throw std::runtime_error(
"Invalid function '" + name +
"(" + args[0] +
"," + args[1] +
")'.");
1020 std::string argsStr =
"(";
1021 for(
size_t i = 0; i < args.size(); i++) {
1022 if(i > 0) argsStr +=
",";
1026 throw std::runtime_error(
"Invalid function/arguments for '" + name + argsStr +
"'.");
1029 template <
typename Derived>
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) {
1038 std::string firstStr = str.substr(0, pos);
1039 std::string lastStr = str.substr(pos + 1);
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));
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;
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);
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) {
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;
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;
1082 throw std::runtime_error(
"Invalid numeric range '" + str +
"'.");
1087 template <
typename Derived>
1092 else if(str.size() == 2) {
1093 size_t pos = mOperators2.find(str);
1094 return (pos != std::string::npos && pos % 2 == 0);
1099 template <
typename Derived>
1103 # ifdef EIGENLAB_DEBUG 1104 bool operationPerformed =
false;
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();
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();
1125 # ifdef EIGENLAB_DEBUG 1126 operationPerformed =
true;
1132 # ifdef EIGENLAB_DEBUG 1133 if(operationPerformed) { std::cout <<
"i: "; printChunks(chunks); std::cout << std::endl; }
1138 template <
typename Derived>
1142 # ifdef EIGENLAB_DEBUG 1143 bool operationPerformed =
false;
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;
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();
1159 lhs = chunks.erase(op);
1160 op = (lhs != chunks.end()) ? lhs + 1 : lhs;
1161 rhs = (op != chunks.end()) ? op + 1 : op;
1163 # ifdef EIGENLAB_DEBUG 1164 operationPerformed =
true;
1174 # ifdef EIGENLAB_DEBUG 1175 if(operationPerformed) { std::cout <<
"-: "; printChunks(chunks); std::cout << std::endl; }
1180 template <
typename Derived>
1184 # ifdef EIGENLAB_DEBUG 1185 bool operationPerformed =
false;
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 ==
".^")) {
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]);
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]);
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();
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));
1213 lhs->value.mapLocal();
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));
1221 lhs->value.mapLocal();
1224 throw std::runtime_error(
"Invalid operand dimensions for operation '" + lhs->field + op->field + rhs->field +
"'.");
1226 chunks.erase(op, rhs + 1);
1228 rhs = (op != chunks.end()) ? op + 1 : op;
1230 # ifdef EIGENLAB_DEBUG 1231 operationPerformed =
true;
1241 # ifdef EIGENLAB_DEBUG 1242 if(operationPerformed) { std::cout <<
"^: "; printChunks(chunks); std::cout << std::endl; }
1247 template <
typename Derived>
1251 # ifdef EIGENLAB_DEBUG 1252 bool operationPerformed =
false;
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 ==
"./")) {
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]);
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]);
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);
1274 lhs->value.local().array() /= rhs->value.matrix()(0, 0);
1276 if(op->field ==
"*" || op->field ==
".*")
1277 lhs->value.local() = lhs->value.matrix().array() * rhs->value.matrix()(0, 0);
1279 lhs->value.local() = lhs->value.matrix().array() / rhs->value.matrix()(0, 0);
1280 lhs->value.mapLocal();
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;
1288 lhs->value.local() = Derived::Constant(rhs->value.matrix().rows(), rhs->value.matrix().cols(), temp).array() / rhs->value.matrix().array();
1289 lhs->value.mapLocal();
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();
1296 lhs->value.local().array() /= rhs->value.matrix().array();
1298 if(op->field ==
".*")
1299 lhs->value.local() = lhs->value.matrix().array() * rhs->value.matrix().array();
1301 lhs->value.local() = lhs->value.matrix().array() / rhs->value.matrix().array();
1302 lhs->value.mapLocal();
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();
1310 lhs->value.local() = lhs->value.matrix() * rhs->value.matrix();
1311 lhs->value.mapLocal();
1315 throw std::runtime_error(
"Invalid operand dimensions for operation '" + lhs->field + op->field + rhs->field +
"'.");
1317 chunks.erase(op, rhs + 1);
1319 rhs = (op != chunks.end()) ? op + 1 : op;
1321 # ifdef EIGENLAB_DEBUG 1322 operationPerformed =
true;
1332 # ifdef EIGENLAB_DEBUG 1333 if(operationPerformed) { std::cout <<
"*: "; printChunks(chunks); std::cout << std::endl; }
1338 template <
typename Derived>
1342 # ifdef EIGENLAB_DEBUG 1343 bool operationPerformed =
false;
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 ==
".-")) {
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]);
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]);
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);
1365 lhs->value.local().array() -= rhs->value.matrix()(0, 0);
1367 if(op->field ==
"+" || op->field ==
".+")
1368 lhs->value.local() = lhs->value.matrix().array() + rhs->value.matrix()(0, 0);
1370 lhs->value.local() = lhs->value.matrix().array() - rhs->value.matrix()(0, 0);
1371 lhs->value.mapLocal();
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;
1379 lhs->value.local() = Derived::Constant(rhs->value.matrix().rows(), rhs->value.matrix().cols(), temp).array() - rhs->value.matrix().array();
1380 lhs->value.mapLocal();
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();
1387 lhs->value.local().array() -= rhs->value.matrix().array();
1389 if(op->field ==
"+" || op->field ==
".+")
1390 lhs->value.local() = lhs->value.matrix().array() + rhs->value.matrix().array();
1392 lhs->value.local() = lhs->value.matrix().array() - rhs->value.matrix().array();
1393 lhs->value.mapLocal();
1397 throw std::runtime_error(
"Invalid operand dimensions for operation '" + lhs->field + op->field + rhs->field +
"'.");
1399 chunks.erase(op, rhs + 1);
1401 rhs = (op != chunks.end()) ? op + 1 : op;
1403 # ifdef EIGENLAB_DEBUG 1404 operationPerformed =
true;
1414 # ifdef EIGENLAB_DEBUG 1415 if(operationPerformed) { std::cout <<
"+: "; printChunks(chunks); std::cout << std::endl; }
1420 template <
typename Derived>
1424 # ifdef EIGENLAB_DEBUG 1425 bool operationPerformed =
false;
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();) {
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]);
1437 if(lhs->type ==
VALUE) {
1438 lhs->value.local() = rhs->value.matrix();
1439 lhs->value.mapLocal();
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();
1447 mVariables[lhs->field].local() = rhs->value.matrix();
1448 mVariables[lhs->field].mapLocal();
1451 lhs->value.matrix().block(lhs->row0, lhs->col0, lhs->rows, lhs->cols) = rhs->value.matrix();
1454 mVariables[lhs->field].local() = rhs->value.matrix();
1455 mVariables[lhs->field].mapLocal();
1458 rhs = chunks.erase(op, rhs + 1);
1459 op = (rhs != chunks.begin()) ? rhs - 1 : rhs;
1460 if (op != chunks.begin()) lhs = op - 1;
1462 # ifdef EIGENLAB_DEBUG 1463 operationPerformed =
true;
1473 # ifdef EIGENLAB_DEBUG 1474 if(operationPerformed) { std::cout <<
"=: "; printChunks(chunks); std::cout << std::endl; }
1480 # ifdef EIGENLAB_DEBUG 1481 template <
typename Derived>
1485 for(
typename ChunkArray::iterator it = chunks.begin(); it != chunks.end(); it++) {
1488 std::cout << textRepresentation(it->value, maxRows, maxCols, precision);
1490 std::cout <<
"(" << it->row0 <<
":" << it->row0 + it->rows - 1 <<
"," << it->col0 <<
":" << it->col0 + it->cols - 1 <<
")";
1493 std::cout << it->field;
1495 std::cout <<
"(" << it->row0 <<
":" << it->row0 + it->rows - 1 <<
"," << it->col0 <<
":" << it->col0 + it->cols - 1 <<
")";
1498 std::cout << it->field;
1501 std::cout <<
"f()=" << it->field;
1508 template <
typename Derived>
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;
1515 template <
typename Derived>
1518 if(val.
matrix().size() == 1)
1519 return numberToString<typename Derived::Scalar>(val.
matrix()(0, 0), precision);
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);
1528 str += (val.
matrix().cols() > maxCols ?
"...]" :
"]");
1530 str += (val.
matrix().rows() > maxRows ?
"...]" :
"]");
1534 # endif // #ifdef EIGENLAB_DEBUG 1535 #endif // #ifdef DEBUG 1537 template <
typename Derived>
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);
1548 template <
typename Derived>
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)));
1559 args.push_back(std::string(i0, str.end()));
1563 template <
typename Derived>
1564 template <
typename T>
1567 std::istringstream iss(str);
1574 return (!iss.fail() && !iss.bad() && iss.eof());
1577 template <
typename Derived>
1578 template<
typename T>
1581 std::istringstream iss(str);
1584 if(iss.fail() || iss.bad() || !iss.eof())
1585 throw std::runtime_error(
"Failed to convert " + str +
" to a number.");
1589 template <
typename Derived>
1590 template<
typename T>
1593 std::ostringstream oss;
1595 oss << std::setprecision(precision) << num;
1602 template <
typename Derived>
1605 typename Derived::Scalar & ,
1609 Derived & , std::true_type)
1615 Derived resultMatrix;
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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);
1667 if(resultMatrix.isApprox(resultValue.
matrix())) std::cout <<
"OK" << std::endl;
1668 else { std::cout <<
"FAIL" << std::endl; ++numFails; }
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);
1678 if(resultMatrix.isApprox(resultValue.
matrix())) std::cout <<
"OK" << std::endl;
1679 else { std::cout <<
"FAIL" << std::endl; ++numFails; }
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; }
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; }
1694 template <
typename Derived>
1697 typename Derived::Scalar & ,
1701 Derived & , std::false_type)
1706 template <
typename Derived>
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;
1715 size_t numFails = 0;
1717 Derived resultMatrix;
1719 typename Derived::Scalar
s = 2;
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);
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);
1737 std::cout <<
"Testing basic operations..." << std::endl << std::endl;
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
1884 std::cout <<
"Test matrix coefficient-wise power a .^ b: ";
1885 resultValue =
eval(
"abs(a) .^ b");
1887 for(
Index i = 0; i < a34.size(); ++i)
1888 resultMatrix(i) = pow(std::abs(a34(i)), b34(i));
1894 if(resultMatrix.isApprox(resultValue.
matrix())) std::cout <<
"OK" << std::endl;
1895 else { std::cout <<
"FAIL" << std::endl; ++numFails; }
1897 std::cout <<
"Test matrix/scalar coefficient-wise power a ^ s: ";
1898 resultValue =
eval(
"abs(a) ^ s");
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; }
1905 std::cout <<
"Test scalar/matrix coefficient-wise power s ^ b: ";
1906 resultValue =
eval(
"s ^ b");
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; }
1913 std::cout <<
"Test matrix/scalar coefficient-wise power a .^ s: ";
1914 resultValue =
eval(
"abs(a) .^ s");
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; }
1921 std::cout <<
"Test scalar/matrix coefficient-wise power s .^ b: ";
1922 resultValue =
eval(
"s .^ b");
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; }
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; }
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; }
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; }
1948 std::cout << std::endl <<
"Testing coefficient and submatrix block access..." << std::endl << std::endl;
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
1999 std::cout <<
"Test subvector segment access using ':' identifier v(:): ";
2000 resultValue =
eval(
"v(:)");
2002 if(resultMatrix.isApprox(resultValue.
matrix())) std::cout <<
"OK" << std::endl;
2003 else { std::cout <<
"FAIL" << std::endl; ++numFails; }
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; }
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; }
2018 std::cout << std::endl <<
"Testing vector/matrix expressions..." << std::endl << std::endl;
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; }
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; }
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; }
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; }
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; }
2058 std::cout << std::endl <<
"Testing coefficient-wise functions..." << std::endl << std::endl;
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
2123 std::cout << std::endl <<
"Testing matrix reduction functions..." << std::endl << std::endl;
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
2207 std::cout << std::endl <<
"Testing matrix functions..." << std::endl << std::endl;
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; }
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; }
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; }
2229 std::cout << std::endl <<
"Testing matrix initializers..." << std::endl << std::endl;
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; }
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; }
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; }
2251 std::cout <<
"Test zeros(5.2): ";
2252 resultValue =
eval(
"zeros(5.2)");
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;
2260 std::cout <<
"Test eye(-3): ";
2261 resultValue =
eval(
"eye(-3)");
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;
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; }
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; }
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; }
2287 std::cout << std::endl <<
"Testing variable assignment..." << std::endl << std::endl;
2290 std::cout <<
"Test assigning to a variable with the same dimensions a = b: ";
2292 if(a34.isApprox(b34)) std::cout <<
"OK" << std::endl;
2293 else { std::cout <<
"FAIL" << std::endl; ++numFails; }
2295 std::cout <<
"Test assigning to a variable with different dimensions 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);
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; }
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; }
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; }
2317 std::cout <<
"Test bad function call: ";
2318 resultValue =
eval(
"foobar(-3)");
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;
2325 std::cout << std::endl <<
"Testing fp parsing..." << std::endl << std::endl;
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; }
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; }
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; }
2346 std::cout << std::endl;
2348 std::cout <<
"Successfully completed unit test for EigenLab with no failures." << std::endl;
2350 std::cout <<
"Completed unit test for EigenLab with " << numFails <<
" failures (see above)." << std::endl;
2351 std::cout << std::endl;
2354 #endif // #ifdef DEBUG 2358 #endif // #ifndef EigenLab_H void evalMatrixExpression(const std::string &str, Value< Derived > &mat)
std::map< std::string, Value< Derived > > ValueMap
void evalFunction(const std::string &name, std::vector< std::string > &args, Value< Derived > &result)
void setShared(const Value &val)
void setShared(const Derived &mat)
const Eigen::Map< Derived > & matrix() const
std::vector< std::string > mFunctions
static std::vector< std::string > split(const std::string &str, const char delimeter)
void setShared(const typename Derived::Scalar *data, size_t rows=1, size_t cols=1)
typename std::is_same< bool, decltype(test< T >(0))>::type type
Value< Eigen::MatrixXd > ValueXd
void clearCachedExpressions()
void setLocal(const Eigen::MatrixBase< Derived > &mat)
void setLocal(const typename Derived::Scalar *data, size_t rows=1, size_t cols=1)
Value< Derived > & var(const std::string &name)
Chunk(const std::string &str="", int t=-1, const Value< Derived > &val=Value< Derived >())
Parser< Eigen::MatrixXd > ParserXd
bool hasVar(const std::string &name)
void evalNumericRange(const std::string &str, Value< Derived > &mat)
void evalIndexRange(const std::string &str, int *first, int *last, int numIndices)
bool evalFunction_2_lt(const std::string &name, Value< Derived > &arg0, Value< Derived > &arg1, Value< Derived > &result, std::false_type)
Eigen::Map< Derived > & matrix()
void evalPowers(ChunkArray &chunks)
std::vector< std::string > splitArguments(const std::string &str, const char delimeter) const
static std::string numberToString(T num, int precision=0)
Value< Eigen::MatrixXf > ValueXf
void clearVar(const std::string &name)
void setCacheExpressions(bool b)
Value(const typename Derived::Scalar s)
bool evalFunction_1_lt(const std::string &name, Value< Derived > &arg, Value< Derived > &result, std::false_type)
void splitEquationIntoChunks(const std::string &expression, ChunkArray &chunks, std::string &code)
Parser< Eigen::MatrixXf > ParserXf
Value< Eigen::MatrixXi > ValueXi
Value & operator=(const typename Derived::Scalar s)
std::map< std::string, ChunkArray > mCachedChunkedExpressions
void evalAddition(ChunkArray &chunks)
const Derived & local() const
static T stringToNumber(const std::string &str)
Parser< Eigen::MatrixXi > ParserXi
static std::string trim(const std::string &str)
bool mCacheChunkedExpressions
bool isOperator(const char c) const
void setLocal(const Value &val)
void evalIndices(ChunkArray &chunks)
bool cacheExpressions() const
void setLocal(const typename Derived::Scalar s)
Value(const Derived &mat)
bool isVariable(const std::string &name) const
void evalAssignment(ChunkArray &chunks)
void evalNegations(ChunkArray &chunks)
std::string::const_iterator findClosingBracket(const std::string &str, const std::string::const_iterator openingBracket, const char closingBracket) const
void evalMultiplication(ChunkArray &chunks)
Value(const typename Derived::Scalar *data, size_t rows=1, size_t cols=1)
bool isFunction(const std::string &str) const
Eigen::Map< Derived > mShared
std::vector< Chunk > ChunkArray
Value< Derived > eval(const std::string &expression)
static bool isNumber(const std::string &str, T *num=0)