　　一年多前做曲线拟合，当时需要用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 = a₀ + a₁*x + a₂*x²。换成线性方程组就是

[1, X₁, X₁²] [a₀] [y₁]

[1, X₂, X₂²] *[a₁] = [y₂]

…… [a₂]

[1, X_n, X_n²] [y_n]

　　在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使用。

发表于 2017-03-09 21:23 Siriuslzx 阅读(3251) 评论(0) 编辑收藏举报

VC++调用R语言

一、用R作曲线拟合

二、C++调用R来实现

三、gcc编译成库给vc++用

一、配置

二、编码

三、使用

三、gcc编译成库给vc++用　