我正在尝试从rtype
对象转换为boost
中的ublas
。
使用我从Rcpp dev列表中找到的一些关于ublas的代码,我能够返回包装为rtype
的ublas
向量。
例如
// Converts from ublas to rtype
template <typename T>
Rcpp::Vector< Rcpp::traits::r_sexptype_traits<T>::rtype >
ublas2rcpp( const boost::numeric::ublas::vector<T>& x ){
return Rcpp::Vector< Rcpp::traits::r_sexptype_traits<T>::rtype >(
x.begin(), x.end()
) ;
}
为了模拟行为导入行为,我目前在相应长度的ublas
向量上使用一个循环,并将从rtype
到它的所有内容都分配给它
从循环切换会提高性能吗?
// My attempt to convert from rtype to ublas
// so R can find the libraries
//[[Rcpp::depends(BH)]]
//[[Rcpp::plugins("cpp11")]]
#include <Rcpp.h>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <boost/numeric/odeint.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
using namespace std;
using namespace Rcpp;
using namespace boost::numeric::ublas;
typedef boost::numeric::ublas::vector< double > vector_type;
typedef boost::numeric::ublas::matrix< double > matrix_type;
template <typename T> /// this need to be fixed up, hopefully working for now though
Rcpp::Vector< Rcpp::traits::r_sexptype_traits<T>::rtype >
ublas2rcpp( const boost::numeric::ublas::vector<T>& x ){
return Rcpp::Vector< Rcpp::traits::r_sexptype_traits<T>::rtype >(
x.begin(), x.end()
) ;
}
// [[Rcpp::export]]
NumericVector main(NumericVector x1)
{
int L =x1.length();
vector_type x(L , 0 ); // initialize the vector to all zero
for(int i=0;i<L;i++)
{
x(i) = x1(i);
}
return(ublas2rcpp(x));
}
首先,很抱歉延迟。希望我所做的一切能弥补。
话虽如此,让我们这样做吧!
引言
问题有两个方面:
- 从R到C++的转换(
Rcpp::as<T>(obj)
) - 从C++转换为R(
Rcpp::wrap(obj)
)
幸运的是,有一个很棒的Rcpp小插曲,名为Extending Rcpp,用于处理自定义对象。遗憾的是,与其他文档相比,小插曲的清晰度还有很多不足之处。(我可能会尝试做一个PR来改善它。)
所以,我将试着用一些评论带你走过这些步骤。请注意,所使用的方法是通过模板和部分专业化,最终会产生一些不错的automagic。
解释
第1阶段-远期申报
在第一阶段,在使用Rcpp.h
之前,我们必须声明我们想要使用的功能。为此,我们将加载一个不同的头文件,并向Rcpp::traits
命名空间添加一些定义。
原则上,当我们开始编写文件时,我们必须加载的第一个标头是RcppCommon.h
,而不是通常的Rcpp.h
!!如果我们不在Rcpp.h
调用之前放置转发声明,我们将无法适当地注册我们的扩展。
然后,我们必须为sourceCpp()
提供不同的插件标记,以便在代码编译期间设置适当的标志。在插件之后,我们将包含我们想要使用的实际boost头。最后,我们必须在Rcpp::traits
命名空间中添加两个特殊的Rcpp函数声明Rcpp::as<T>(obj)
和Rcpp::wrap(obj)
。要启用多个类型,我们必须创建一个Exporter
类,而不是更直接地调用template <> ClassName as( SEXP )
。
#include <RcppCommon.h>
// Flags for C++ compiler
// [[Rcpp::depends(BH)]]
// [[Rcpp::plugins("cpp11")]]
// Third party library includes that provide the template class of ublas
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/matrix.hpp>
// Provide Forward Declarations
namespace Rcpp {
namespace traits{
// Setup non-intrusive extension via template specialization for
// 'ublas' class boost::numeric::ublas
// Support for wrap
template <typename T> SEXP wrap(const boost::numeric::ublas::vector<T> & obj);
// Support for as<T>
template <typename T> class Exporter< boost::numeric::ublas::vector<T> >;
}
}
第2阶段-包括Rcpp.h
仅仅有一个阶段来申报进口订单可能看起来很无聊,但如果在前向申报之前包括Rcpp.h
,那么Rcpp::traits
就不会更新,我们就进入了深渊。
因此:
// >> Place <Rcpp.h> AFTER the forward declaration!!!! <<
#include <Rcpp.h>
// >> Place Definitions of Forward Declarations AFTER <Rcpp.h>!!!! <<
第3阶段-实施扩展
现在,我们必须切实执行前瞻性声明。特别地,唯一稍微有问题的实现是as<>
,因为wrap()
是直接的。
wrap()
要实现wrap()
,我们必须调用Rcpp中名为Rcpp::traits::r_sexptype_traits<T>::rtype
的内置类型转换索引。由此,我们能够获得包含RTYPE
的int
,然后构建Rcpp::Vector
。对于矩阵的构造,同样的想法也是成立的。
as()
对于as<>()
,我们需要考虑将要传入的模板。此外,我们在Exporter
类定义的正下方设置了一个typedef
,以便轻松定义要在get()
方法中使用的OUT
对象。除此之外,我们使用相同的技巧从C++T
类型来回移动到R
类型。
为了完成as<>
,或者从R到C++的直接端口,我不得不做一些肮脏的事情:我复制了向量内容。控制该输出的代码在Exporter
类的get()
中给出。您可能希望花一些时间来研究使用指针来更改分配。我不太熟悉ublas
,所以我没有看到一种简单的方法来解决指针传递问题
// Define template specializations for as<> and wrap
namespace Rcpp {
namespace traits{
// Defined wrap case
template <typename T> SEXP wrap(const boost::numeric::ublas::vector<T> & obj){
const int RTYPE = Rcpp::traits::r_sexptype_traits<T>::rtype ;
return Rcpp::Vector< RTYPE >(obj.begin(), obj.end());
};
// Defined as< > case
template<typename T>
class Exporter< boost::numeric::ublas::vector<T> > {
typedef typename boost::numeric::ublas::vector<T> OUT ;
// Convert the type to a valid rtype.
const static int RTYPE = Rcpp::traits::r_sexptype_traits< T >::rtype ;
Rcpp::Vector<RTYPE> vec;
public:
Exporter(SEXP x) : vec(x) {
if (TYPEOF(x) != RTYPE)
throw std::invalid_argument("Wrong R type for mapped 1D array");
}
OUT get() {
// Need to figure out a way to perhaps do a pointer pass?
OUT x(vec.size());
std::copy(vec.begin(), vec.end(), x.begin()); // have to copy data
return x;
}
} ;
}
}
第4阶段-测试
好吧,让我们看看我们所做的工作是否得到了回报(扰流板成功了!扰流器)。要检查,我们应该看两个不同的领域:
- 在函数和中跟踪诊断
- 自动魔术测试
下面给出了这两种情况。注意,我选择将ublas
设置缩短为:
// Here we define a shortcut to the boost ublas class to enable multiple ublas types via a template.
// ublas::vector<T> => ublas::vector<double>, ... , ublas::vector<int>
namespace ublas = ::boost::numeric::ublas;
跟踪诊断
// [[Rcpp::export]]
void containment_test(Rcpp::NumericVector x1) {
Rcpp::Rcout << "Converting from Rcpp::NumericVector to ublas::vector<double>" << std::endl;
ublas::vector<double> x = Rcpp::as< ublas::vector<double> >(x1); // initialize the vector to all zero
Rcpp::Rcout << "Running output test with ublas::vector<double>" << std::endl;
for (unsigned i = 0; i < x.size (); ++ i)
Rcpp::Rcout << x(i) << std::endl;
Rcpp::Rcout << "Converting from ublas::vector<double> to Rcpp::NumericVector" << std::endl;
Rcpp::NumericVector test = Rcpp::wrap(x);
Rcpp::Rcout << "Running output test with Rcpp::NumericVector" << std::endl;
for (unsigned i = 0; i < test.size (); ++ i)
Rcpp::Rcout << test(i) << std::endl;
}
测试呼叫:
containment_test(c(1,2,3,4))
结果:
Converting from Rcpp::NumericVector to ublas::vector<double>
Running output test with ublas::vector<double>
1
2
3
4
Converting from ublas::vector<double> to Rcpp::NumericVector
Running output test with Rcpp::NumericVector
1
2
3
4
此测试按预期进行。进入下一次测试!
Automagic测试
// [[Rcpp::export]]
ublas::vector<double> automagic_ublas_rcpp(ublas::vector<double> x1) {
return x1;
}
测试呼叫:
automagic_ublas_rcpp(c(1,2,3.2,1.2))
结果:
[1] 1.0 2.0 3.2 1.2
成功!
第5阶段-所有Cntrl+C和<kbd]Cntrl>+Vp>以下是stage给出的上述代码块的组合。如果将其复制并粘贴到.cpp
文件中,则的所有内容都应正常工作。如果没有,请告诉我。
// -------------- Stage 1: Forward Declarations with `RcppCommon.h`
#include <RcppCommon.h>
// Flags for C++ compiler
// [[Rcpp::depends(BH)]]
// [[Rcpp::plugins("cpp11")]]
// Third party library includes that provide the template class of ublas
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/matrix.hpp>
// Here we use ublas_vec to enable multiple ublas types via a template.
// ublas::vector<T> => ublas::vector<double>, ... , ublas::vector<int>
namespace ublas = ::boost::numeric::ublas;
// Provide Forward Declarations
namespace Rcpp {
namespace traits{
// Setup non-intrusive extension via template specialization for
// 'ublas' class boost::numeric::ublas
// Support for wrap
template <typename T> SEXP wrap(const boost::numeric::ublas::vector<T> & obj);
// Support for as<T>
template <typename T> class Exporter< boost::numeric::ublas::vector<T> >;
}
}
// -------------- Stage 2: Including Rcpp.h
// >> Place <Rcpp.h> AFTER the forward declaration!!!! <<
#include <Rcpp.h>
// >> Place Definitions of Forward Declarations AFTER <Rcpp.h>!!!! <<
// -------------- Stage 3: Implementation of Declarations
// Define template specializations for as<> and wrap
namespace Rcpp {
namespace traits{
// Defined wrap case
template <typename T> SEXP wrap(const boost::numeric::ublas::vector<T> & obj){
const int RTYPE = Rcpp::traits::r_sexptype_traits<T>::rtype ;
return Rcpp::Vector< RTYPE >(obj.begin(), obj.end());
};
// Defined as< > case
template<typename T>
class Exporter< boost::numeric::ublas::vector<T> > {
typedef typename boost::numeric::ublas::vector<T> OUT ;
// Convert the type to a valid rtype.
const static int RTYPE = ::Rcpp::traits::r_sexptype_traits< T >::rtype ;
Rcpp::Vector<RTYPE> vec;
public:
Exporter(SEXP x) : vec(x) {
if (TYPEOF(x) != RTYPE)
throw std::invalid_argument("Wrong R type for mapped 1D array");
}
OUT get() {
// Need to figure out a way to perhaps do a pointer pass?
OUT x(vec.size());
std::copy(vec.begin(), vec.end(), x.begin()); // have to copy data
return x;
}
} ;
}
}
// -------------- Stage 4: Tests
// [[Rcpp::export]]
ublas::vector<double> automagic_ublas_rcpp(ublas::vector<double> x1) {
return x1;
}
// [[Rcpp::export]]
void containment_test(Rcpp::NumericVector x1) {
Rcpp::Rcout << "Converting from Rcpp::NumericVector to ublas::vector<double>" << std::endl;
ublas::vector<double> x = Rcpp::as< ublas::vector<double> >(x1); // initialize the vector to all zero
Rcpp::Rcout << "Running output test with ublas::vector<double>" << std::endl;
for (unsigned i = 0; i < x.size (); ++ i)
Rcpp::Rcout << x(i) << std::endl;
Rcpp::Rcout << "Converting from ublas::vector<double> to Rcpp::NumericVector" << std::endl;
Rcpp::NumericVector test = Rcpp::wrap(x);
Rcpp::Rcout << "Running output test with Rcpp::NumericVector" << std::endl;
for (unsigned i = 0; i < test.size (); ++ i)
Rcpp::Rcout << test(i) << std::endl;
}
结束语
Whew。。。太多了。希望以上提供了足够的理由,因为我相信您可能希望将1D
向量扩展到ublas::matrix
,依此类推。此外,等待应该是值得的,因为您现在拥有Rcpp
的自动转换魔力,因此无需在return()
语句中调用ublas2rcpp()
。事实上,您可以将函数的返回类型简化为ublas::vector<double>
!
在其他新闻中,我认为可能阻止您创建上面的备用模板函数(例如rcpp2ublas
)的原因是无法推断Rcpp::Vector
是什么C++T
类型。就我个人而言,在Rcpp GitHub Repo中挖掘后,我不确定是否存在这样的转换索引,也不确定由于Shield
的使用,是否应该使用它。深入挖掘转化率是另一天的冒险。
// -------------- Stage 1: Forward Declarations with `RcppCommon.h`
#include <RcppCommon.h>
// Flags for C++ compiler
// [[Rcpp::depends(BH)]]
// [[Rcpp::plugins("cpp11")]]
// Third party library includes that provide the template class of ublas
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/matrix.hpp>
// Here we use ublas_vec to enable multiple ublas types via a template.
// ublas::vector<T> => ublas::vector<double>, ... , ublas::vector<int>
namespace ublas = ::boost::numeric::ublas;
// Provide Forward Declarations
namespace Rcpp {
namespace traits{
// Setup non-intrusive extension via template specialization for
// 'ublas' class boost::numeric::ublas
// Support for wrap
template <typename T> SEXP wrap(const boost::numeric::ublas::vector<T> & obj);
// Support for as<T>
template <typename T> class Exporter< boost::numeric::ublas::vector<T> >;
}
}
// -------------- Stage 2: Including Rcpp.h
// >> Place <Rcpp.h> AFTER the forward declaration!!!! <<
#include <Rcpp.h>
// >> Place Definitions of Forward Declarations AFTER <Rcpp.h>!!!! <<
// -------------- Stage 3: Implementation of Declarations
// Define template specializations for as<> and wrap
namespace Rcpp {
namespace traits{
// Defined wrap case
template <typename T> SEXP wrap(const boost::numeric::ublas::vector<T> & obj){
const int RTYPE = Rcpp::traits::r_sexptype_traits<T>::rtype ;
return Rcpp::Vector< RTYPE >(obj.begin(), obj.end());
};
// Defined as< > case
template<typename T>
class Exporter< boost::numeric::ublas::vector<T> > {
typedef typename boost::numeric::ublas::vector<T> OUT ;
// Convert the type to a valid rtype.
const static int RTYPE = ::Rcpp::traits::r_sexptype_traits< T >::rtype ;
Rcpp::Vector<RTYPE> vec;
public:
Exporter(SEXP x) : vec(x) {
if (TYPEOF(x) != RTYPE)
throw std::invalid_argument("Wrong R type for mapped 1D array");
}
OUT get() {
// Need to figure out a way to perhaps do a pointer pass?
OUT x(vec.size());
std::copy(vec.begin(), vec.end(), x.begin()); // have to copy data
return x;
}
} ;
}
}
// -------------- Stage 4: Tests
// [[Rcpp::export]]
ublas::vector<double> automagic_ublas_rcpp(ublas::vector<double> x1) {
return x1;
}
// [[Rcpp::export]]
void containment_test(Rcpp::NumericVector x1) {
Rcpp::Rcout << "Converting from Rcpp::NumericVector to ublas::vector<double>" << std::endl;
ublas::vector<double> x = Rcpp::as< ublas::vector<double> >(x1); // initialize the vector to all zero
Rcpp::Rcout << "Running output test with ublas::vector<double>" << std::endl;
for (unsigned i = 0; i < x.size (); ++ i)
Rcpp::Rcout << x(i) << std::endl;
Rcpp::Rcout << "Converting from ublas::vector<double> to Rcpp::NumericVector" << std::endl;
Rcpp::NumericVector test = Rcpp::wrap(x);
Rcpp::Rcout << "Running output test with Rcpp::NumericVector" << std::endl;
for (unsigned i = 0; i < test.size (); ++ i)
Rcpp::Rcout << test(i) << std::endl;
}
1D
向量扩展到ublas::matrix
,依此类推。此外,等待应该是值得的,因为您现在拥有Rcpp
的自动转换魔力,因此无需在return()
语句中调用ublas2rcpp()
。事实上,您可以将函数的返回类型简化为ublas::vector<double>
!rcpp2ublas
)的原因是无法推断Rcpp::Vector
是什么C++T
类型。就我个人而言,在Rcpp GitHub Repo中挖掘后,我不确定是否存在这样的转换索引,也不确定由于Shield
的使用,是否应该使用它。深入挖掘转化率是另一天的冒险。