Help testing release candidate for std::complex

Developers
54

@Bob_Carpenter @ChrisChiasson I see this thread has been inactive for a while, so don’t know if this is something you are still actively working on, or whether discussion has moved to another thread?

In any case the Eigen-AD repository is now publicly available linked to from here.

I spent a bit of time trying to reproduce the issues you have been having. Initially I thought that swapping Eigen 3.3.3 for Eigen-AD resolves these issues, but now it seems that my little toy example works with both versions of Eigen. I.e. I can multiply a matrix of double variables with a vector of std::complex<stan::math::var> variables.

I would be interested to know whether there are still unresolved issues and whether you are planning to merge this feature branch into the development branch anytime soon?

Here is a complete example that works on my machine (with Eigen 3.3.3),

// Note that prototype for supporting complex arithmetic of Stan types is currently in the following branch
// https://github.com/stan-dev/math/tree/feature/0123-complex-var
#include <stan/math.hpp>

template<typename BinaryOp>
struct Eigen::ScalarBinaryOpTraits<std::complex<stan::math::var>, double, BinaryOp> {
  typedef std::complex<stan::math::var> ReturnType;
};

template<typename BinaryOp>
struct Eigen::ScalarBinaryOpTraits<double, std::complex<stan::math::var>, BinaryOp> {
  typedef std::complex<stan::math::var> ReturnType;
};

typedef Eigen::Matrix<stan::math::var, Eigen::Dynamic, Eigen::Dynamic> Stan_MT;
typedef Eigen::Matrix<std::complex<stan::math::var>, Eigen::Dynamic, Eigen::Dynamic> Stan_CMT;
typedef Eigen::Matrix<std::complex<stan::math::var>, Eigen::Dynamic, 1> Stan_CVT;

using namespace std;
using namespace Eigen;

int main(){
  try{
    Stan_MT A = Stan_MT::Random(6,6);
    cout << "Here is a random 6x6 matrix, A:" << endl << A << endl << endl;

    EigenSolver<Stan_MT> es(A);
    cout << "The eigenvalues of A are:" << endl << es.eigenvalues() << endl << endl;

    Stan_CMT V = es.eigenvectors(); // this line causes a compiler error
    cout << "The matrix of eigenvectors, V, is:" << endl << es.eigenvectors() << endl << endl;

    complex<stan::math::var> lambda = es.eigenvalues()[0];
    cout << "Consider the first eigenvalue, lambda = " << lambda << endl;

    Stan_CVT v = es.eigenvectors().col(0);
    cout << "If v is the corresponding eigenvector, then lambda * v = " << endl << lambda * v << endl;

    Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> A_double = Eigen::MatrixXd::Random(6,6);
    Stan_CVT x1_cvt = A_double * v;

  } catch (const std::exception& e){
    std::cout << " a standard exception was caught, with message '"
              << e.what() << "'\n";
  }

  return 0;
}