标签: armadillo

在我的代码中,我有循环,我构建和确定线性系统,并尝试解决它:

#pragma omp parallel for
for (int i = 0; i < n[0]+1; i++) {
    for (int j = 0; j < n[1]+1; j++) {
        for (int k = 0; k < n[2]+1; k++) {
            arma::mat A(max_points, 2);
            arma::mat y(max_points, 1);
            // initialize A and y

            arma::vec solution = solve(A,y);
        }
    }
}

有时,程序会随机挂起,或者解决方案向量中的结果是NaN.如果我这样做:

arma::vec solution;
#pragma omp critical 
{
    solution = solve(weights*A,weights*y);
}

然后这些问题似乎不再发生了.

当它挂起时,它会这样做,因为有些线程正在等待OpenMP屏障:

Thread 2 (Thread 0x7fe4325a5700 (LWP 39839)):
#0  0x00007fe44d3c2084 in gomp_team_barrier_wait_end () from /usr/lib64/gcc-4.9.2/lib64/gcc/x86_64-redhat-linux-gnu/4.9.2/libgomp.so.1
#1 …

考虑以下R代码,

## ----------- R version -----------

caller <- function(x=1:3, fun = "identity", ...){

  ## do some other stuff
  ## ...
  ## then call the function
  eval(call(fun, x))

}

fun1 <- function(x, ...){
  x + x
}

fun2 <- function(x, a = 10) a * x

caller(fun = "fun1")
caller(fun = "fun2")

用户可以传递函数名称"fun" caller.我希望用RcppArmadillo对象执行相同的任务(显然是一个更复杂的任务的一部分).该函数将被定义C++,并且用户通过引用其名称在R级别选择它:

caller_cpp(1:3, "fun1_cpp")

要么

caller_cpp(1:3, "fun2_cpp")

等等

这是我对调用函数的天真尝试,甚至无法编译:

## ----------- C++ version -----------

library(Rcpp)
require( RcppArmadillo )    

sourceCpp( code = '

       // [[Rcpp::depends("RcppArmadillo")]] …

我在blitz ++,armadillo,boost :: MultiArray之间用以下代码进行了比较(借用旧帖子)

#include <iostream>
using namespace std;
#include <windows.h>
#define _SCL_SECURE_NO_WARNINGS
#define BOOST_DISABLE_ASSERTS 
#include <boost/multi_array.hpp>
#include <blitz/array.h>
#include <armadillo>

int main(int argc, char* argv[])
{
    const int X_SIZE = 1000;
    const int Y_SIZE = 1000;
    const int ITERATIONS = 100;
    unsigned int startTime = 0;
    unsigned int endTime = 0;

    // Create the boost array


    //------------------Measure boost Loop------------------------------------------
    {
        typedef boost::multi_array<double, 2> ImageArrayType;
        ImageArrayType boostMatrix(boost::extents[X_SIZE][Y_SIZE]);
        startTime = ::GetTickCount();
        for (int i = 0; i < ITERATIONS; …

我的任务是用C++重写R函数来加速while循环.所有R代码都在Rcpp和Armadillo的帮助下重写,除了.Fortran().我尝试首先使用Rinside,它的工作速度非常慢,正如Dirk所指出的那样.(数据通过R - > C++ - > R - > Fortran是昂贵的)

由于我不想用C++重写Fortran代码,反之亦然,通过将C++直接链接到Fortran来加速程序看起来很自然:R - > C++ - > Fortran.

// [[Rcpp::depends(RcppArmadillo)]]

#include <RcppArmadillo.h>
using namespace Rcpp;

extern "C"{
   List f_(int *n,NumericMatrix a, NumericVector c, double* eps);
}

问题是我可以将C++与Fortran集成并将R与C++集成,但我不能让这三个东西一起工作!

我尝试在Linux中编译C++,但它无法找到RcppArmadillo.h并且namespace Rcpp:

 error: RcppArmadillo.h: No such file or directory
 error: 'Rcpp' is not a namespace-name

当我sourceCpp("test.cpp")直接打电话给R时,控制台会显示:

test.o:test.cpp:(.text+0x20b2): undefined reference to `f_'
collect2: ld returned 1 exit status
Error in sourceCpp("test.cpp") : Error occurred building shared library. …

我正在尝试将armadillo库安装到我的linux系统(ubuntu 12.04)上.首先需要BOOST BLAS ATLAS和LAPACK进行安装.有没有办法检查这些库是否已经安装？

仅仅为了完成我的C++/Rcpp编程,我开始实现(样本)标准偏差函数:

#include <Rcpp.h>
#include <vector>
#include <cmath>
#include <numeric>

// [[Rcpp::export]]
double cppSD(Rcpp::NumericVector rinVec)
{
  std::vector<double> inVec(rinVec.begin(),rinVec.end());
  int n = inVec.size();
  double sum = std::accumulate(inVec.begin(), inVec.end(), 0.0);
  double mean = sum / inVec.size();

  for(std::vector<double>::iterator iter = inVec.begin();
      iter != inVec.end(); ++iter){
        double temp;
        temp= (*iter - mean)*(*iter - mean);
        *iter = temp;
      }

  double sd = std::accumulate(inVec.begin(), inVec.end(), 0.0);
  return std::sqrt( sd / (n-1) );
}

我还决定测试stddevArmadillo库中的函数,考虑到它可以在向量上调用:

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;

// [[Rcpp::export]]

double armaSD(arma::colvec …

C中的任何函数是否等同%in%于R中的运算符？考虑R中的以下命令:

which(y %in% x)

我试图在C++中找到相同的东西(特别是在Armadillo中),我找不到任何东西.然后我编写了自己的函数,与上面的R命令相比,这个函数非常慢.

这是我写的:

#include <RcppArmadillo.h>
// [[Rcpp::depends("RcppArmadillo")]]

// [[Rcpp::export]]
arma::uvec myInOperator(arma::vec myBigVec, arma::vec mySmallVec ){
 arma::uvec rslt = find(myBigVec == mySmallVec[0]);
 for (int i = 1; i < mySmallVec.size(); i++){
   arma::uvec rslt_tmp = find(myBigVec == mySmallVec[i]);
   rslt = arma::unique(join_cols( rslt, rslt_tmp ));
 }
 return rslt;
}

现在,在上面的代码中采购后,我们有:

x <- 1:4
y <- 1:10
res <- benchmark(myInOperator(y, x), which(y %in% x), columns = c("test",
      "replications", "elapsed", "relative", "user.self", "sys.self"), 
       order = "relative")

以下是结果:

                 test replications elapsed …

当给出相同的输入时,R sum()函数和RcppArmadillo accu()函数的结果存在细微差别.例如,以下代码:

R:

vec <- runif(100, 0, 0.00001)
accu(vec)
sum(vec)

C++:

// [[Rcpp::depends("RcppArmadillo")]]
// [[Rcpp::export]]
double accu(arma::vec& obj)
{
    return arma::accu(obj);
}

给出结果:

0.00047941851844312633(C++)

0.00047941851844312628(R)

根据http://keisan.casio.com/calculator的说法,真正的答案是:

4.79418518443126270948E -4-

这些小的差异在我的算法中加起来并且显着影响它的执行方式.有没有办法在C++中更准确地总结向量？或者至少在不必调用R代码的情况下获得与R相同的结果？

我是新手使用犰狳,尽管尝试/搜索相当多,但无法得到以下内容.

我需要执行两个巨大的(动态)数组(不是向量).我决定使用犰狳.我理解如何使用向量初始化arma :: mat,但是我可以使用数组吗？我不明白,因为我在文档中没有看到任何提及.我试图避免出于内部设计原因使用向量.我尝试使用样本数组手动初始化每个元素(作为一个愚蠢但起点).像下面的代码不起作用.

using namespace std;
using namespace arma;   

mat A(SIZE, 1), B(SIZE, 1);

for(int i = 0; i < SIZE; i++)
{
    A << v[i] << endr;
    B << c[i] << endr;
}

cout << "A: " << endl;
A.print();
cout << "B: " << endl;
B.print();

对于输入数组v = {1,2,0,-1,.9}和c = {0,.5,1,-2,-5}.输出将是:

A:
        0
B:
  -5.0000

这是可以理解的.任何使用数组初始化arma :: mat或arma :: colvector的工作？提前致谢!

我希望我的重写问题现在符合Stackoverflow的标准.请考虑以下示例.我正在编写一个Log-Likelihood函数,其中计算cdf over vectors是最耗时的部分.示例1使用R::pnorm,示例2近似于正常cdf erfc.正如您所看到的结果非常相似,ercf版本更快一点.

实际上(在MLE中)然而事实证明,ercf并不精确,这使得算法可以进入inf区域,除非准确地设置约束.我的问题:

1)我错过了什么吗？是否有必要实现一些错误处理(对于erfc)？

2)您是否有任何其他建议来加快代码或替代方案？考虑并行化for循环是否有回报？

require(Rcpp)
require(RcppArmadillo)
require(microbenchmark)

#Example 1 : standard R::pnorm
src1 <- '
NumericVector ppnorm(const arma::vec& x,const arma::vec& mu,const     arma::vec& sigma, int lt, int lg) {
int n = x.size();
arma::vec res(n);
for (int i=0; i<n; i++) {
   res(i) = R::pnorm(x(i),mu(i),sigma(i),lt,lg);
}
return wrap(res);
}
'

#Example 2: approximation with ercf
src2 <- '
NumericVector ppnorm(const arma::vec& x,const arma::vec& mu,const    arma::vec& sigma, int lt, int lg) {
int n = x.size();
arma::vec …