一年多前做曲线拟合，当时需要用C++调用R语言来完成。
一、用R作曲线拟合
先看一段用R语言作拟合的示例：
x <- runif(100,min=0,max=100) #创建100个随机数 y <- x*x+runif(x,-10,10)*x+10*rnorm(x) #创建y向量 plot(x,y) #绘制散点图 matr <- data.frame(X=x, X2=x*x) #建立解析矩阵 fm <- lsfit(matr ,y) #最小二乘法解析 a <- coef(fm) #从解析结果提取系数 f <- function(x) a[1]+ a["X"]*x+a["X2"]*x*x #建立拟合函数 curve(f,col="red",add=TRUE) #画曲线
运行结果为：
这里创建的数据接近二次函数，因此拟合的目标是y = a0 + a1*x + a2*x2。换成线性方程组就是
[1, X1, X12] [a0] [y1]
[1, X2, X22] *[a1] = [y2]
…… [a2]
[1, Xn, Xn2] [yn]
在R语言中，如果是求解该线性方程组，可以使用solve(Matr,y)；如果是用最小二乘法拟合，就要用函数lsfit()，它返回的结果再提取系数就可以得到系数的行列式。
二、C++调用R来实现
如果用c++调用R来做，就需要借助工具。我使用的是R语言下的两个包，Rcpp与RInside，作用分别是R调用C++与C++调用R，其中RInside依赖于Rcpp。RInside提供的接口如下：
int parseEval(const std::string &line, SEXP &ans); // parse line, return in ans; error code rc void parseEvalQ(const std::string &line); // parse line, no return (throws on error) void parseEvalQNT(const std::string &line); // parse line, no return (no throw) Proxy parseEval(const std::string &line); // parse line, return SEXP (throws on error) Proxy parseEvalNT(const std::string &line); // parse line, return SEXP (no throw) template <typename T> void assign(const T& object, const std::string& nam) { global_env_m->assign( nam, object ) ; } …… Rcpp::Environment::Binding operator[]( const std::string& name );
其中parseEval函数是用来执行语句的，而assign与[]是用来赋值的。
所以使用RInside，以上过程就是这样：
1 RInside RI; 2 RI["M"] = Rcpp::NumericMatrix(/*传入矩阵*/); 3 RI["y"] = Rcpp::NumericVector(/*传入数组*/); 4 std::vector<double> b2 = RI.parseEval("coef(lsfit(M,y))");
最初是直接使用RInside的，结果编译报错，说是不支持VC++，没办法，只好先用gcc编译成库，再给VC++使用。
三、gcc编译成库给vc++用
一、配置
首先下载R，Rtools（使用里面的gcc），Rcpp，RInside，CodeBlocks（不带编译器），在CodeBlocks的compiler settings中如下设置：
PS：那时候Rtools里的gcc还是4.6.3，不支持C++11，为了使用C++11，费了好大功夫，最后失败。现在gcc版本跟上来了。
然后在CodeBlocks中新建一个dll工程：
设置工程属性：
以下3张图分别设置链接库，头文件目录，库目录，对应着gcc的-l、-I、-L命令。
PS：user32其实是不必要的。
二、编码
这里只展示两个简单的函数，最小二乘拟合与解线性方程组，非常简单。
头文件：
1 #ifndef RINSIDE_4_VS_H 2 #define RINSIDE_4_VS_H 3 4 #include <boost/container/vector.hpp> 5 #include <boost/multi_array.hpp> 6 typedef boost::container::vector<double> vector_t; 7 typedef boost::multi_array<double, 2> matrix_t; 8 9 /** 10 * 11 * author Li zhongxian 12 * version 1.0 13 * date 2016.2 14 * 15 */ 16 #ifdef BUILD_DLL 17 #define DLL_EXPORT __declspec(dllexport) 18 #else 19 #define DLL_EXPORT __declspec(dllimport) 20 #endif 21 22 #ifdef __cplusplus 23 extern "C" 24 { 25 #endif 26 /** rief 对于线性方程组M·b=y，用最小二乘法拟合 27 * 28 * param M 解析矩阵，大小为(系数-1)*案例 29 * param y 指标向量，大小为案例数目 30 * param b 系数向量，大小为系数数目 31 * eturn 异常信息 32 */ 33 DLL_EXPORT const char* lsfit(const matrix_t &M, const vector_t &y, vector_t &b); 34 DLL_EXPORT const char* lsfit_c(double M[], size_t rows, size_t cols, double y[], double b[]); 35 /** rief 对于线性方程组M·b=y，求解b 36 * 37 * param M 样本矩阵，大小为系数*案例 38 * param y 指标向量，大小为案例数目 39 * param b 系数向量，大小为系数数目 40 * eturn 异常信息 41 */ 42 DLL_EXPORT const char* solve(const matrix_t &M, const vector_t &y, vector_t &b); 43 DLL_EXPORT const char* solve_c(double M[], size_t rows, size_t cols, double y[], double b[]); 44 45 #ifdef __cplusplus 46 } 47 #endif 48 49 #endif //H
源文件：
1 #include "RInside4vs.h" 2 #include <RInside.h> 3 //#include <fstream> 4 5 inline Rcpp::NumericMatrix createMatrix(const matrix_t &Mat) 6 { 7 int rows = Mat.shape()[0]; 8 int cols = Mat.shape()[1]; 9 Rcpp::NumericMatrix M(cols,rows,Mat.data()); 10 return M; 11 } 12 13 inline Rcpp::NumericMatrix createMatrix_c(double data[], size_t rows, size_t cols) 14 { 15 Rcpp::NumericMatrix M(cols,rows,data); 16 return M; 17 } 18 19 /// 所有函数共用一个RInside对象 20 static RInside RI; 21 22 DLL_EXPORT const char* lsfit(const matrix_t &M, const vector_t &y, vector_t &b) 23 { 24 25 // std::ostringstream sout; 26 // std::ofstream fs("Rlog.txt"); 27 // std::streambuf *x = std::cout.rdbuf(fs.rdbuf()); 28 // char msg[200]=""; 29 try{ 30 RI["M"] = createMatrix(M); 31 // RI.parseEvalQ("cat('Showing M '); print(M)"); 32 RI["y"] = y; 33 // RI.parseEvalQ("cat('Showing y '); print(y)"); 34 35 std::vector<double> b2 = RI.parseEval("coef(lsfit(M,y))"); 36 std::copy(b2.begin(), b2.end(), b.begin()); 37 38 }catch(std::exception &e){ 39 // std::cerr << "failed in " << __FUNCTION__ << std::endl; 40 // return strcat(strcpy(msg,sout.str().c_str()), e.what()); 41 return e.what(); 42 } 43 /// 清理R空间 44 RI.parseEvalQ("rm(list = ls())"); 45 // std::cerr.rdbuf(x); 46 return ""; 47 } 48 49 DLL_EXPORT const char* lsfit_c(double M[], size_t rows, size_t cols, double y[], double b[]) 50 { 51 52 try{ 53 RI["M"] = createMatrix_c(M,rows,cols); 54 // RI.parseEvalQ("cat('Showing M '); print(M)"); 55 RI["y"] = Rcpp::NumericVector(y, y+rows); 56 // RI.parseEvalQ("cat('Showing y '); print(y)"); 57 58 std::vector<double> b2 = RI.parseEval("coef(lsfit(M,y))"); 59 std::copy(b2.begin(), b2.end(), b); 60 61 }catch(std::exception &e){ 62 return e.what(); 63 } 64 /// 清理R空间 65 RI.parseEvalQ("rm(list = ls())"); 66 return ""; 67 } 68 69 DLL_EXPORT const char* solve(const matrix_t &M, const vector_t &y, vector_t &b) 70 { 71 72 // std::cout << "len(y): " << y.size() << std::endl; 73 // std::cout << "len(M): " << M.size() << std::endl; 74 // std::cout << "rows: " << rows << std::endl; 75 // std::cout << "cols: " << cols << std::endl; 76 77 try{ 78 RI["M"] = createMatrix(M); 79 // RI.parseEvalQ("cat('Showing M '); print(M)"); 80 RI["y"] = y; 81 // RI.parseEvalQ("cat('Showing y '); print(y)"); 82 83 std::vector<double> b2 = RI.parseEval("solve(M,y)"); 84 std::copy(b2.begin(), b2.end(), b.begin()); 85 86 }catch(std::exception &e){ 87 return e.what(); 88 } 89 /// 清理R空间 90 RI.parseEvalQ("rm(list = ls())"); 91 return ""; 92 } 93 94 DLL_EXPORT const char* solve_c(double M[], size_t rows, size_t cols, double y[], double b[]) 95 { 96 97 try{ 98 RI["M"] = createMatrix_c(M,rows,cols); 99 // RI.parseEvalQ("cat('Showing M '); print(M)"); 100 RI["y"] = Rcpp::NumericVector(y, y+rows); 101 // RI.parseEvalQ("cat('Showing y '); print(y)"); 102 103 std::vector<double> b2 = RI.parseEval("solve(M,y)"); 104 std::copy(b2.begin(), b2.end(), b); 105 106 }catch(std::exception &e){ 107 return e.what(); 108 } 109 /// 清理R空间 110 RI.parseEvalQ("rm(list = ls())"); 111 return ""; 112 }
三、使用
编译完成后，会生成libRInside4vs.a、libRInside4vs.def、RInside4vs.dll，只有dll是需要用到的。.a是gcc使用的.lib，VS不能用。打开cmd，切换到RInside4vs.dll所在的目录，执行下面两条语句：
pexports RInside4vs.dll -o > RInside4vs.def
lib /def:RInside4vs.def
其中，pexports是gcc的扩展工具的命令，需要下载mingw-utils-0.3。lib是VS的命令。它们的配置这里就不细讲了。
执行完成后，即生成VS可用的RInside4vs.lib，这里的lib是引导文件，并不是静态库，它配合RInside4vs.dll使用。