c++ Rcpp & main/optim.c -通过`nmmin()`传递变量后如何翻译`void *ex`?

7ivaypg9  于 2023-07-01  发布在  其他
关注(0)|答案(1)|浏览(122)

**主要问题:**如何在Rcpp中格式化一个复杂的结构体,使其能被optim.c接受?当返回到定义的optimfn函数时,如何修改此结构体?
**背景:**我正在Rcpp中编写一个EM算法,将用于beta混合模型。我试图通过R API(nmmin)使用一个简单的nelder mead数值优化方法,但是我很难通过变量列表。我对C++/Rcpp非常陌生。
信息:

eout计数:

  • zprobs:矩阵(double[45203 x5])。(列和行的数量将变化)。
  • parmlist:一个包含avec、mvec和svec的列表(每个列表有5个double)。(长度会有所不同)。
  • xi:长度为45203的双精度数列表。(长度会有所不同)。
  • denom:长度为45203的双精度数列表。(长度会有所不同)。

为了将我的变量传递给nmmin,我尝试将它转换为结构体格式,首先将NumericMatrix和NumericVector转换为double *var,(但我后来了解到这些不是真实的(编辑:这些不包含向量)),然后用这些项创建一个具有我定义的结构的对象。当前的方法能够编译,但我不确定*ex进入enlikeB函数时是什么样子。我还尝试将Rcpp::List项转换为std::vectors,以防这些项可能通过。

我现在在哪里:

#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);`
 }
z5btuh9x

z5btuh9x1#

解决了!
我定义了一个类来代替typedef struct

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) {}
};

在将数据转换为正确的格式后,我能够使void *ex对象等于:

Ething et(xid, xil, zrows, zcols, zdata, denomd, denoml);

这允许您在创建对象时超越设置内存/大小,并以易于解释的方式传递格式化的对象。完整的代码将在几个月内在github上的mgaynor1/nQuack下提供。

相关问题