R与Cpp混合编程二
Last update: 2021-06-09
Rcpp中自动支持缺失值
find_nonfinite可以查找NA，NaN, infinite.
arma::field<uvec> obsindex(arma::mat Xmis){
int p = Xmis.n_cols, j;
field<uvec> indexfield(p);
for(j= 0; j< p; ++j){
indexfield(j) = find_nonfinite(Xmis.col(j));
}
return(indexfield);
}List返回值的类型的重复申明
- 为什么cpp中定义的函数返回List后，List的成分为矩阵时，无法做index_max()计算和返回该矩阵到另一个List的成分中？？
返回的List的成分为矩阵时，需要重新申明其为矩阵。
List aList = getList(xx);
mat mat1 = aList["mat"]; // 注意这里的申明mat不能省略，否则无法强制转换。
注意事项
若A是一个矩阵，则A.col(k)可以直接支持diagmat(),repmat()函数中使用，表示它返回一个列向量。
若C是一个cube，则C.slice(k)是一个矩阵，支持矩阵的所有运算，比如.t(),.i()以及A + b*C.slice(k);而 C.col(k)不是一个矩阵而是一个subcube，不支持直接的矩阵运算，只能通过赋值给一个矩阵tmpMat后转化为一个矩阵后，再进行相应的矩阵运算。为了更加高效的使用cube，请注意cube的存储格式，将需要参与矩阵乘法的subcube储存到第三个维度上。
下面展示field
和cube运算的例子，其中C是一个qqK的cube，B是一个K元的field ，它的每一个成分都是一个矩阵。
C.slice(k).i() * B(r).col(j) + B(r).col(l)- 最大值最小值溢出，C++中最小值。非常接近0的溢出为0。比如\(4^{-200}*A %( 2* 4^{-200}*B)\)，对于这种就需要先把\(4^{-200}\)约掉后再计算两个矩阵的逐元除法。
2)整除错误：C++中整数除整数等于相除之后再向下取整。如下示例代码：
U0 = (bx - 1/n *A * tA * a) / as_scalar(1- 1/n * a.t()*tA * a);
## 这里1/n=0 而不等于真正的小数。应该修改成：
U0 = (bx - 1.0/n *A * tA * a) / as_scalar(1- 1.0/n * a.t()*tA * a);我知道这样C++有这个语法规则，但是写代码的时候还是容易疏忽，而且还找不到错误的点在哪里。
- 公式中的加号运算错误的写成减法运算，导致出错。这是我容易出错的地方。 例如：
X_tk = X - repmat(U0.t() + Mu0.row(k)* W0.t(), n, 1);
#错误写成：
X_tk = X - repmat(U0.t() - Mu0.row(k)* W0.t(), n, 1);- 避免log(a)当a=0报错的情况，可以使用
log(a + (a==0)) # a=0时,值会等于0如何传递非结构化数据
R中list类型常常用于存储非结构化数据，如何将非结构化数据传入到C++中进行计算，这是一个非常重要的问题。 下面这段使用List的代码非常经典，值得学习。
#include "RcppArmadillo.h"
// [[Rcpp::depends(RcppArmadillo)]]
#include<ctime>
#define INT_MIN (-INT_MAX - 1)
using namespace Rcpp;
using namespace arma;
using namespace std;
//[[Rcpp::export]]
List calList(const Rcpp::List& alist, const Rcpp::List blist){
int r, L1 = alist.length();
cout<<"L1="<<L1<<endl;
field<mat> Xf(1,2), Wf(1,2);
for(r=0; r < L1; ++r){
mat Xtmp = alist[r];
Xf(0,r) = Xtmp;
mat Wtmp = blist[r];
Wf(0,r) = Wtmp;
}
List resList = List::create(
Rcpp::Named("Xf") = Xf,
Rcpp::Named("Wf") = Wf
);
return resList;
}类型转换
field<mat> Rf1 = ICM_fit["Rf"];
field<ivec> yf(r_max);
uvec y1_u = index_max(Rf1(r),1);
yf(r) = conv_to< ivec >::from(y1_u) + 1;逐行逐列做运算
使用.each_col()和.each_row()属性的用法，能够避免使用循环。
mat X(6, 5, fill::ones);
vec v = linspace<vec>(10,15,6);
X.each_col() += v; // in-place addition of v to each column vector of X
mat Y = X.each_col() + v; // generate Y by adding v to each column vector of X
// subtract v from columns 0 through to 3 in X
X.cols(0,3).each_col() -= v;
uvec indices(2);
indices(0) = 2;
indices(1) = 4;
X.each_col(indices) = v; // copy v to columns 2 and 4 in XList套List
Armadillo中List套List的经典代码:
//[[Rcpp::export]]
List solution_pathCpp(const vec& Y, const mat& X,const vec& lambda_set, const double& alpha0, const double& r2,
mat& beta_int, const int& maxIter=100, const double& epsObj=1e-4, const bool& verbose=true){
int i, nLam = lambda_set.n_elem, p = X.n_cols;
mat betaMat(p ,nLam, fill::zeros);
List resList(nLam);
List tmpList = lmPenAlphaCpp(Y, X, alpha0, lambda_set(0), r2, beta_int, maxIter,
epsObj, false, verbose);
mat betai = tmpList["beta0"];
resList[0] = tmpList;
betaMat.col(0) = betai;
if(nLam > 1){
for(i=1; i< nLam; ++i){
tmpList = lmPenAlphaCpp(Y, X, alpha0, lambda_set(0), r2, betai, maxIter,
epsObj, true, verbose);
vec tmpvec = tmpList["beta0"];
betaMat.col(i) = tmpvec;
resList[i] = tmpList;
}
}
List outList = List::create(
Rcpp::Named("betaMat") = betaMat,
Rcpp::Named("lambda_set") = lambda_set,
Rcpp::Named("fitList") = resList
);
return(outList);
}非规则向量元素访问
Armadillo中的向量原来还可以这样访问元素 gamma(S)：
uvec S = find(abs(gamma) <= r2);
cout<<"S="<<S.t()<<endl;
cout<<gamma(S).t()<<endl;
double obj = pow(norm(Y- X* (gamma+ alpha)), 2.0) / n;
if(S.n_elem >0){
pen2 = lambda2* accu(gamma(S) % gamma(S)) * pow(S.n_elem, -nu);
}else{
pen2 = 1e20;
}