Main question: How to format a complex struct in Rcpp that will be accepted by optim.c? and how will this struct be modified when returned to the defined optimfn function?
Background: I'm writing an EM algorithm in Rcpp that will be used in a beta mixture model. I am trying to use a simple nelder mead numerical optimization approach through the R API (nmmin), however I am having a hard time passing the variable list through. I'm very very new to C++/Rcpp.
Information:
eout countains:
To pass my variables through to nmmin, I've tried to convert it to a struct format by first converting the NumericMatrix and NumericVector into double *var, (but I then learned these are not real (edit: that these do not hold vectors)), followed by creating an object with my defined struct with these items. The current method is able to compile, but I am not sure what *ex looks like when it enters the enlikeB function. I also tried to convert the Rcpp::List items to std::vectors, in case these were possible to pass through.
Where I am currently:
#include <RcppCommon.h>
typedef struct{
double *xi;
int xil;
int zrows;
int zcols;
double *zdata;
double *denom;
int denoml;
}Ething;
#include <Rcpp.h>
//HELPER FUNCTIONS HERE:
//
// Defining the optimfn function:
typedef double optimfn(int n, double *par, void *ex);
double elnlikeB(int n, double *par, void *ex) {
Rcpp::Rcout << "entering elnlikeB\n";
Rcpp::Rcout << "The value of ex: " << ex << "\n";
Ething *et = (Ething *) ex;
Rcpp::Rcout << "The value of et : " << et << "\n";
Rcpp::Rcout << "making z\n";
int zrows = et ->zrows;
int zcols = et ->zcols;
double *z = et ->zdata;
Rcpp::NumericMatrix Zprobsmat;
int count = 0;
for(int x = 0; x < zcols; x++) {
for(int y = 0; y < zrows; y++) {
Zprobsmat(y, x) = z[count]; //row, column
count = count +1;
}
}
Rcpp::Rcout << "The value of zprob : " << Zprobsmat << "\n";
Rcpp::Rcout << "making xi\n";
int xil = et -> xil;
double *xid = et -> xi;
Rcpp::NumericVector xi;
for(int nn = 0; nn < xil; nn++){
xi[nn] = xid[nn];
}
Rcpp::Rcout << "The value of xi : " << xi << "\n";
Rcpp::Rcout << "making denom\n";
int denoml = et ->denoml;
double *denomd = et -> denom;
Rcpp::NumericVector denom;
for(int nn = 0; nn < denoml; nn++){
denom[nn] = denomd[nn];
}
Rcpp::Rcout << "The value of denom: " << denom << "\n";
Rcpp::Rcout << "making myguess\n";
Rcpp::NumericVector myguess;
for(int nn = 0; nn < n; nn++){
myguess[nn] = par[nn];
}
Rcpp::Rcout << "The value of myguess: " << myguess << "\n";
//This section works is the vectors and matrix are correct//
Rcpp::Rcout << "calculating\n";
int nnew = xi.size();
int nmixt = Zprobsmat.ncol();
int gsize = myguess.size();
Rcpp::NumericVector avec(nmixt);
avec[Rcpp::Range(0, (nmixt - 2))] = invmlogitc(myguess[Rcpp::Range(0, (nmixt - 2))]);
avec[(nmixt-1)] = (1 - rcpp_sum(avec));
Rcpp::NumericVector mvec(nmixt);
Rcpp::NumericVector svec(nmixt);
mvec = mutransform(myguess[Rcpp::Range((nmixt-1), (((nmixt-1)*2)))]);
svec = stransform(myguess[Rcpp::Range((((nmixt-1)*2)+1), (gsize-1))]);
Rcpp::NumericMatrix loglikemat( nnew, nmixt );
Rcpp::NumericMatrix lnfr( nnew, nmixt );
double sum = 0;
for(int i = 0; i < nmixt; i++){
double a = log(avec[i]);
double m = mvec[i];
double s = svec[i];
Rcpp::NumericVector ab = alphabetacalc(m, s);
for( int z = 0; z < nnew; z++){
lnfr(z,i) = (R::dbeta( xi[z], ab[0], ab[1], true));
double left = (Zprobsmat(z,i)/denom[z]);
double right = (a + lnfr(z,i));
loglikemat(z, i) = (left*right);
sum += loglikemat(z, i);
}
}
sum = (-sum);
Rprintf("sum: %d \n", sum);
return(sum);
}
optimfn elnlikeB;
//More helpers
extern "C" {
void nmmin(int n, double *xin, double *x, double *Fmin, optimfn fn,
int *fail, double abstol, double intol, void *ex,
double alpha, double beta, double gamma, int trace,
int *fncount, int maxit);
}
// [[Rcpp::export]]
Rcpp::List optim_test(Rcpp::List eout){
Rcpp::Rcout << "starting\n";
Rcpp::Rcout << "making avec, mvec, svec\n";
Rcpp::List parmlist = eout["parm.list"];
Rcpp::NumericVector avec = parmlist["avec"];
Rcpp::NumericVector mvec = parmlist["mvec"];
Rcpp::NumericVector svec = parmlist["svec"];
Rcpp::Rcout << "making zprob and xi \n";
Rcpp::NumericMatrix Zprobsmat = eout["zprob"];
Rcpp::NumericVector xi = eout["xi"];
int xil = xi.size();
Rcpp::NumericVector denom = eout["denom"];
int denoml = denom.size();
double nmixt = avec.size();
Rcpp::NumericVector mas = mlogitc(avec, nmixt);
Rcpp::Rcout << "making guess\n";
int gl = ((nmixt-1) + nmixt + nmixt);
Rcpp::NumericVector guess(gl);
guess[Rcpp::Range(0, (nmixt-2))] = mas;
guess[Rcpp::Range((nmixt-1), ((nmixt-1)*2))] = min(mvec);
guess[Rcpp::Range((((nmixt-1)*2)+1), (gl-1))] = sin(svec);
Rcpp::Rcout << "making zdata\n";
int tl = (xil*nmixt);
Rcpp::NumericVector zdata(tl);
int count = 0;
for(int x = 0; x < nmixt; x++) {
for(int y = 0; y < xil; y++) {
zdata[count] = Zprobsmat(y, x); //row, column
count = count +1;
}
}
int zrows = xil;
int zcols = nmixt;
Rcpp::Rcout << "making et\n";
Rcpp::List et = Rcpp::List::create(Rcpp::Named("xi") = xi, Rcpp::Named("xil") = xil,
Rcpp::Named("zrows") = zrows,
Rcpp::Named("zcols") = zcols,
Rcpp::Named("zdata") = zdata, Rcpp::Named("denom") = denom,
Rcpp::Named("denoml") =denoml);
Rcpp::Rcout << "making guess vec\n";
double vec[gl];
for(int nn = 0; nn < gl; nn++){
vec[nn] = guess[nn];
}
double opar[gl];
double Fmin = 0.0;
int fail = 0;
double abstol = 1.0e-8;
double intol = 1.0e-8;
double alpha = 1.0;
double beta = 0.5;
double gamma = 2.0;
int trace =0;
int fncount = 0;
int maxit = 500;
Rcpp::Rcout << "Attempting nmmin\n";
nmmin(gl, vec, opar, &Fmin,
elnlikeB, &fail, abstol, intol, &et, alpha, beta,
gamma, trace, &fncount, maxit);
Rcpp::List res = Rcpp::List::create(Rcpp::_["Fmin"] = Fmin, Rcpp::_["fail"]=fail);
return(res);`
}
Solved!
Instead of typedef struct, I defined a class:
class Ething {
public:
double *xi;
int xil;
int zrows;
int zcols;
double *zdata;
double *denom;
int denoml;
Ething(double *xi, int xil, int zrows, int zcols, double *zdata, double *denom, int denoml) : xi(xi), xil(xil), zrows(zrows), zcols(zcols), zdata(zdata), denom(denom), denoml(denoml) {}
};
After converting my data into the correct format, I was able to make the void *ex object equal to:
Ething et(xid, xil, zrows, zcols, zdata, denomd, denoml);
This allows you to surpass setting the memory/size when creating the object and passes the object as formatted in an easy-to-interpret manner. The full code will be available on github under mgaynor1/nQuack in a few months.