trying out RcppNumeric (C++ script)

numeric.cpp

/ [[Rcpp::depends(RcppEigen)]]                                                                  
// [[Rcpp::depends(RcppNumerical)]]                                                              
                                                                                                 
#include <RcppNumerical.h>                                                                       
using namespace Numer;                                                                           
                                                                                                 
typedef Eigen::Map<Eigen::MatrixXd> MapMat;                                                      
typedef Eigen::Map<Eigen::VectorXd> MapVec;                                                      
                                                                                                 
class FooFun: public MFuncGrad                                                                   
{                                                                                                
private:                                                                                         
  const MapMat pcoefs;                                                                           
  int i;                                                                                         
public:                                                                                          
  FooFun(const MapMat coefs_, int i_) : pcoefs(coefs_), i(i_) {}                                 
                                                                                                 
  double f_grad(Constvec& beta, Refvec grad)                                                     
  {                                                                                              
                                                                                                 
    const double f = pcoefs(0,i)*pow(beta[0],4.0) + pcoefs(1,i)*pow(beta[0],2.0) + pcoefs(2,i)*b\
eta[0] + pcoefs(3,i)*pow(beta[1],2.0) + pcoefs(4,i);                                             
                                                                                                 
    grad[0] = 4*pcoefs(0,i)*pow(beta[0],3.0) + 2*pcoefs(1,i)*beta[0] + pcoefs(2,i);              
    grad[1] = 2*pcoefs(3,i)*beta[1];                                                             
                                                                                                 
    return f;                                                                                    
  }                                                                                              
}; 

// [[Rcpp::export]]                                                                              
Rcpp::NumericMatrix foo_fun(Rcpp::NumericMatrix coefs)                                           
{                                                                                                
  const MapMat mcoefs = Rcpp::as<MapMat>(coefs);                                                 
  int n = coefs.ncol();                                                                          
  Eigen::MatrixXd betas(2, n);                                                                   
  Eigen::VectorXd beta(2);                                                                       
                                                                                                 
  double fopt;                                                                                   
  for (int i = 0; i < n; i++) {                                                                  
    FooFun foo(mcoefs, i);                                                                       
    beta.setZero();                                                                              
    int status = optim_lbfgs(foo, beta, fopt);                                                   
    betas.col(i) = beta;                                                                         
  }                                                                                              
                                                                                                 
  return Rcpp::wrap(betas);                                                                      
}