Ok, so I’ve tried to just write my function (which is a bit more complex that what I wrote before, however it should only return a simple matrix, and it’s arguments are simple matrices and a real complex number). Here is what I wrote in Stan:
matrix foo(real sReal, real sImag, matrix P, matrix lambda, real Ua, real Ub);
and here is what I wrote in a separate cpp file:
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> foo(const double sReal,const double sImag,
const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> P,
const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> lambda,
const double Ua,
const double Ub,
std::ostream* pstream__) {
std::complex<double> s;
s.real(sReal);
s.imag(sImag);
Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> tmp1 =P.array()*(lambda.array()/(lambda.array()+s));
Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> Id=tmp1;
Id.setIdentity();
Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> tmp2=(Id-tmp1).inverse();
Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> tmp3=tmp2.diagonal().asDiagonal().inverse();
std::complex<double> LTuni = (std::exp(-s*Ua)-std::exp(-s*Ub))/(s*(Ub-Ua));
return((tmp1*tmp2*tmp3*LTuni).real());
}
I then wrote a stupid model with a transformed parameters block that just gets the output of foo for a given set of arguments, and it seems to work - i.e., the results match the real results of this function. Does that mean that I am good to go with implementing the other parts of the function (which do not include complex matrix inversion or complex numbers) in Stan and the model itself? Or am I missing something? I am not sure how exactly does this work with the fact that std::complex should not yet work in Stan…
One issue that I’ve noticed is that I can call this CPP function from the model/transformed parameter/generated quantities block, but if I’m trying to call it from another Stan function - it doesn’t work. For example, I created another Stan function that simply calls foo:
matrix test(real sReal, real sImag, matrix P, matrix lambda, real Ua, real Ub){
matrix[cols(P),cols(P)] a;
a=foo(sReal,sImag,P,lambda,Ua,Ub);
return(a);
}
I get the following error message:
file125c5eca5307.o:file125c5eca5307.cpp:(.text$_ZN60model125c42415a5b_c54d1236bbbfb99ce518497ac68d1123_namespace4testIddddddEEN5Eigen6MatrixIN5boost4math5tools12promote_argsIT_T0_T1_T2_NS6_IT3_T4_ffffE4typeEfE4typeELin1ELin1ELi0ELin1ELin1EEERKS7_RKS8_RKNS2_IS9_Lin1ELin1ELi0ELin1ELin1EEERKNS2_ISA_Lin1ELin1ELi0ELin1ELin1EEERKSB_RKSC_PSo[_ZN60model125c42415a5b_c54d1236bbbfb99ce518497ac68d1123_namespace4testIddddddEEN5Eigen6MatrixIN5boost4math5tools12promote_argsIT_T0_T1_T2_NS6_IT3_T4_ffffE4typeEfE4typeELin1ELin1ELi0ELin1ELin1EEERKS7_RKS8_RKNS2_IS9_Lin1ELin1ELi0ELin1ELin1EEERKNS2_ISA_Lin1ELin1ELi0ELin1ELin1EEERKSB_RKSC_PSo]+0x2d6): undefined reference to `Eigen::Matrix<boost::math::tools::promote_args<double, double, double, doub
What am I missing here?
BW,
Isaac.