一年多前做曲线拟合,当时需要用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     /** \brief 对于线性方程组M·b=y,用最小二乘法拟合
27     *
28     * \param M 解析矩阵,大小为(系数-1)*案例
29     * \param y 指标向量,大小为案例数目
30     * \param b 系数向量,大小为系数数目
31     * \return 异常信息
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     /** \brief 对于线性方程组M·b=y,求解b
36     *
37     * \param M 样本矩阵,大小为系数*案例
38     * \param y 指标向量,大小为案例数目
39     * \param b 系数向量,大小为系数数目
40     * \return 异常信息
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\n'); print(M)");
 32         RI["y"] = y;
 33     //    RI.parseEvalQ("cat('Showing y\n'); 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\n'); print(M)");
 55         RI["y"] = Rcpp::NumericVector(y, y+rows);
 56     //    RI.parseEvalQ("cat('Showing y\n'); 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\n'); print(M)");
 80         RI["y"] = y;
 81     //    RI.parseEvalQ("cat('Showing y\n'); 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\n'); print(M)");
100         RI["y"] = Rcpp::NumericVector(y, y+rows);
101     //    RI.parseEvalQ("cat('Showing y\n'); 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使用。