Test eigen api + possible dynamic lazy wrapper for OpenSCAD

eigentest.cc

# include <Eigen/Core>
# include <iostream>
#include <memory>
/*
g++ -I/usr/include/eigen3 eigentest.cc -o test

*/

// Base virtual class to allow dynamic cast of EigenWrapper<T> even if the various Ts are unrelated
class EigenWrapperBase {
 public:
  virtual ~EigenWrapperBase() {}
};
// Wraps a type T so we can dynamic_cast it.
template <class T>
class EigenWrapper : public EigenWrapperBase {
 public:
  T value;
  EigenWrapper(T&& v) : value(std::move(v)) {}
};

class LazyEigen {
	typedef std::shared_ptr<EigenWrapperBase> EigenPtr;

	EigenPtr intermediatePtr;
	std::function<EigenPtr()> toMatrix;

	template <class T>
	LazyEigen(EigenPtr &&intermediateP, std::function<EigenPtr()> &&toMat)
		//: intermediatePtr(std::move(wrap(std::move(intermediateValue)))),
		: intermediatePtr(std::move(intermediateP)),
			toMatrix(std::move(toMat)) {}

	template <class T, class M>
	LazyEigen(T &&intermediateValue, std::function<M()> &&toMat)
		: LazyEigen(
			std::move(EigenWrapper<T>(std::move(intermediateValue))),
			toMatrix(std::move([toMat=std::move(toMat)]() { return wrap(toMat()); }))) {}
public:

	// template <class T>
	// LazyEigen(T &&intermediateValue, std::function<Matrix3d()> &&toMat);
	// template <class T>
	// LazyEigen(T &&intermediateValue, std::function<Matrix4d()> &&toMat);
	// template <class T>
	// LazyEigen(T &&intermediateValue, std::function<MatrixXd()> &&toMat);

	// template <class T>
	// LazyEigen(T &&intermediateValue, std::function<Vector2d()> &&toMat);
	// template <class T>
	// LazyEigen(T &&intermediateValue, std::function<Vector3d()> &&toMat);
	// template <class T>
	// LazyEigen(T &&intermediateValue, std::function<Vector4d()> &&toMat);
	// template <class T>
	// LazyEigen(T &&intermediateValue, std::function<VectorXd()> &&toMat);

	// template <class T>
	// LazyEigen(EigenPtr &&interm, std::function<EigenPtr()> &&toMat)
	// 	: intermediatePtr(std::move(interm)),
	// 		toMatrix(std::move(toMat)) {}

	template <class T>
	static EigenWrapper<T> wrap(T &&v) {
		return EigenWrapper<T>(std::move(v));
	}
};

#define RECOGNIZE_TYPE_AND_RETURN_UN_OP(operand, type, op) \
  { \
  	if (auto cast = dynamic_pointer_cast<EigenWrapper<type>>(operand->intermediate())) \
		return std::move(LazyEigen( \
			shared_ptr(std::move(op cast)), \
			[=]() { return op cast->matrix(); })); \
  }
#define RECOGNIZE_TYPES_AND_RETURN_BIN_OP(lhs, rhs, lhsType, rhsType, op) \
  { \
  	if (auto lhsCast = dynamic_pointer_cast<EigenWrapper<lhsType>>(lhs->intermediate())) \
  		if (auto rhsCast = dynamic_pointer_cast<EigenWrapper<rhsType>>(rhs->intermediate())) \
  			return std::move(LazyEigen( \
  				shared_ptr(std::move(lhsCast op rhsCast)), \
  				[=]() { return lhsCast->matrix() op rhsCast->matrix(); })); \
  }
#define RECOGNIZE_TYPES_AND_RETURN_UNARY_FUN(operand, type, fn) \
  { \
  	if (auto cast = dynamic_pointer_cast<EigenWrapper<type>>(operand->intermediate())) \
		return std::move(LazyEigen( \
			shared_ptr(std::move(fn(cast))), \
			[=]() { return fn(cast); })); \
  }
#define RECOGNIZE_TYPES_AND_RETURN_BINARY_FUN(lhs, rhs, lhsType, rhsType, fn) \
  { \
  	if (auto lhsCast = dynamic_pointer_cast<EigenWrapper<lhsType>>(lhs->intermediate())) \
  		if (auto rhsCast = dynamic_pointer_cast<EigenWrapper<rhsType>>(rhs->intermediate())) \
  			return std::move(LazyEigen( \
  				shared_ptr(std::move(fn(lhsCast, rhsCast))), \
  				[=]() { return fn(lhsCast->matrix(), rhsCast->matrix()); })); \
  }

/*
typedef Eigen::internal::scalar_sum_op<double, double> ScalarSum;
typedef Eigen::internal::scalar_subtract_op<double, double> ScalarSubtract;
typedef Eigen::internal::scalar_product_op<double, double> ScalarProduct;
typedef Eigen::internal::scalar_divide_op<double, double> ScalarDivide;

/*
 * Some level of dynamic typing to chain efficient Eigen operations.
 * (a + b * (d + e) - f)
 *
 * When lazy eigen values can't be chained in a lazy way, we convert them to matrices.
 *

mult(LazyEigen& lhs, LazyEigen& rhs) {
  // Keep some recognized expression patterns lazy 

  // RECOGNIZE_TYPES_AND_RETURN_BIN_OP(Product<const Matrix4d, const Matrix4d>, const Matrix4d, +);
  // RECOGNIZE_TYPES_AND_RETURN_BIN_OP(CwiseBinaryOp<ScalarSum, const Matrix4d, const Matrix4d>, Vector4d, +);
  // RECOGNIZE_TYPES_AND_RETURN_BIN_OP(Vector4d, CwiseBinaryOp<ScalarSum, const Matrix4d, const Matrix4d>, +);

  // RECOGNIZE_TYPES_AND_RETURN_BIN_OP(const Matrix4d, const Matrix4d, +);
  // RECOGNIZE_TYPES_AND_RETURN_BIN_OP(const Matrix3d, const Matrix3d, +);
  // RECOGNIZE_PRODUCT(Matrix4d, Matrix4d, +);

  // RECOGNIZE_PRODUCT(const Matrix3d, const Matrix3d, +);
  // RECOGNIZE_PRODUCT(const MatrixXd, const MatrixXd, +);
  // RECOGNIZE_PRODUCT(const MatrixXd, const Matrix3Xd, +);
  // RECOGNIZE_PRODUCT(const Matrix3Xd, const Matrix3d, +);

  // RECOGNIZE_CWISE_BIN_OP(ScalarSum, const Matrix4d, const Matrix4d, +);
  // RECOGNIZE_CWISE_BIN_OP(ScalarSubtract, const Matrix4d, const Matrix4d, +);
  // RECOGNIZE_CWISE_BIN_OP(ScalarProduct, const Matrix4d, const Matrix4d, +);
  // RECOGNIZE_CWISE_BIN_OP(ScalarDivide, const Matrix4d, const Matrix4d, +);

  // RECOGNIZE_CWISE_BIN_OP(ScalarSum, const Matrix3d, const Matrix3d, +);
  // RECOGNIZE_CWISE_BIN_OP(ScalarSubtract, const Matrix3d, const Matrix3d, +);
  // RECOGNIZE_CWISE_BIN_OP(ScalarProduct, const Matrix3d, const Matrix3d, +);
  // RECOGNIZE_CWISE_BIN_OP(ScalarDivide, const Matrix3d, const Matrix3d, +);

  // RECOGNIZE_CWISE_BIN_OP(ScalarSum, const MatrixXd, const MatrixXd, +);
  // RECOGNIZE_CWISE_BIN_OP(ScalarSubtract, const MatrixXd, const MatrixXd, +);
  // RECOGNIZE_CWISE_BIN_OP(ScalarProduct, const MatrixXd, const MatrixXd, +);
  // RECOGNIZE_CWISE_BIN_OP(ScalarDivide, const MatrixXd, const MatrixXd, +);

  if (auto lhs1 = dynamic_pointer_cast<CwiseBinaryOp<ScalarSum, const Matrix4d, const Matrix4d>>(expr.intermediate())) {
  	if (auto rhs1 = dynamic_pointer_cast<CwiseBinaryOp<ScalarSum, const Matrix4d, const Matrix4d>>(expr.intermediate())) {
		return LazyEigen(shared_ptr(std::move(lhs1 * rhs1)), [=]() { return lhs.matrix() * rhs.matrix(); });
	}
  }
  auto res = shared_ptr<Matrix4d>(lhs.matrix() * rhs.matrix());
  return LazyEigen(res, [=]() { return res; });
}
*/


using namespace std ;
using namespace Eigen ;

int main () {
	cout << " Eigen version : " << EIGEN_MAJOR_VERSION << " . "
	<< EIGEN_MINOR_VERSION << endl ;

	Matrix3d A ;
	Matrix4d B ;

	// Initialize A
	A << 1.0f , 0.0f , 0.0f ,
		0.0f , 1.0f , 0.0f ,
		0.0f , 0.0f , 1.0f ;

	// Initialize B by accessing individual elements
	for (auto i = 1; i < 4; i++) {
		for (auto j = 1; j < 4; j++) {
			B(j, i ) = 0.0;
		}
	}

	Matrix4d M1 = Matrix4d::Random () ;
	Matrix4d M2 = Matrix4d::Constant (2.2) ;
	// Addition
	// The size and the coefficient - types of the matrices must match
	cout << M1 + M2 << endl ;
	// Matrix multiplication
	// The inner dimensions and the coefficient - types must match
	cout << M1 * M2 << endl ;
	// Scalar multiplication , and subtraction
	// What do you expect the output to be?
	cout << M2 - Matrix4d::Ones () * 2.2 << endl ;

	// Square each element of the matrix
	cout << M1 . array () . square () << endl ;
	// Multiply two matrices element - wise
	cout << M1 . array () * Matrix4d :: Identity () . array () << endl ;
	// All relational operators can be applied element - wise
	// cout << M1 . array () <= M2 . array () << endl << endl ;
	// cout << M1 . array () > M2 . array () << endl ;

	Product<Matrix4d, Matrix4d> prod = M1 * M2;
	typedef Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<double, double>, const Matrix4d, const Matrix4d> MatrixCwiseSum;
	const MatrixCwiseSum expr2 = M2 + M2;
	auto ptr = make_shared<MatrixCwiseSum>(std::move(expr2));
	auto eign = make_shared<LazyEigen>(
			std::move(expr2),
			[=]() { return expr2; });
	cout << "prod: " << prod << "\n";
	cout << "expr2: " << expr2 << "\n";

	return 0;
}