Rcpp Version 0.12.14
Matrix.h
Go to the documentation of this file.
1 // -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
2 //
3 // Matrix.h: Rcpp R/C++ interface class library -- matrices
4 //
5 // Copyright (C) 2010 - 2016 Dirk Eddelbuettel and Romain Francois
6 //
7 // This file is part of Rcpp.
8 //
9 // Rcpp is free software: you can redistribute it and/or modify it
10 // under the terms of the GNU General Public License as published by
11 // the Free Software Foundation, either version 2 of the License, or
12 // (at your option) any later version.
13 //
14 // Rcpp is distributed in the hope that it will be useful, but
15 // WITHOUT ANY WARRANTY; without even the implied warranty of
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 // GNU General Public License for more details.
18 //
19 // You should have received a copy of the GNU General Public License
20 // along with Rcpp. If not, see <http://www.gnu.org/licenses/>.
21 
22 #ifndef Rcpp__vector__Matrix_h
23 #define Rcpp__vector__Matrix_h
24 
25 namespace Rcpp{
26 
27 template <int RTYPE, template <class> class StoragePolicy = PreserveStorage >
28 class Matrix : public Vector<RTYPE, StoragePolicy>, public MatrixBase<RTYPE, true, Matrix<RTYPE,StoragePolicy> > {
29  int nrows ;
30 
31 public:
32  using Vector<RTYPE, StoragePolicy>::size; // disambiguate diamond pattern for g++-6 and later
33 
34  struct r_type : traits::integral_constant<int,RTYPE>{} ;
41 
42  typedef StoragePolicy<Matrix> Storage ;
44  typedef typename VECTOR::iterator iterator ;
47  typedef typename VECTOR::stored_type stored_type ;
48  typedef typename VECTOR::Proxy Proxy ;
49  typedef typename VECTOR::const_Proxy const_Proxy ;
50 
51  Matrix() : VECTOR(Dimension(0, 0)), nrows(0) {}
52 
53  Matrix(SEXP x) : VECTOR( r_cast<RTYPE>( x ) ), nrows( VECTOR::dims()[0] ) {}
54 
55  Matrix( const Dimension& dims) : VECTOR( Rf_allocMatrix( RTYPE, dims[0], dims[1] ) ), nrows(dims[0]) {
56  if( dims.size() != 2 ) throw not_a_matrix();
57  VECTOR::init() ;
58  }
59  Matrix( const int& nrows_, const int& ncols) : VECTOR( Dimension( nrows_, ncols ) ),
60  nrows(nrows_)
61  {}
62 
63  template <typename Iterator>
64  Matrix( const int& nrows_, const int& ncols, Iterator start ) :
65  VECTOR( start, start + (static_cast<R_xlen_t>(nrows_)*ncols) ),
66  nrows(nrows_)
67  {
68  VECTOR::attr( "dim" ) = Dimension( nrows, ncols ) ;
69  }
70 
71  Matrix( const int& n) : VECTOR( Dimension( n, n ) ), nrows(n) {}
72 
73 
74  Matrix( const Matrix& other) : VECTOR( other.get__() ), nrows(other.nrows) {}
75 
76  template <bool NA, typename MAT>
77  Matrix( const MatrixBase<RTYPE,NA,MAT>& other ) : VECTOR( Rf_allocMatrix( RTYPE, other.nrow(), other.ncol() ) ), nrows(other.nrow()) {
78  import_matrix_expression<NA,MAT>( other, nrows, ncol() ) ;
79  }
80 
81  Matrix( const SubMatrix<RTYPE>& ) ;
82 
83  Matrix& operator=(const Matrix& other) {
84  SEXP x = other.get__() ;
85  if( ! ::Rf_isMatrix(x) ) throw not_a_matrix();
86  VECTOR::set__( x ) ;
87  nrows = other.nrows ;
88  return *this ;
89  }
90  Matrix& operator=( const SubMatrix<RTYPE>& ) ;
91 
92  explicit Matrix( const no_init_matrix& obj) {
93  Storage::set__( Rf_allocMatrix( RTYPE, obj.nrow(), obj.ncol() ) );
94  }
95 
96  inline int ncol() const {
97  return VECTOR::dims()[1];
98  }
99  inline int nrow() const {
100  return nrows ;
101  }
102  inline int cols() const {
103  return VECTOR::dims()[1];
104  }
105  inline int rows() const {
106  return nrows ;
107  }
108 
109  inline Row row( int i ){ return Row( *this, i ) ; }
110  inline ConstRow row( int i ) const{ return ConstRow( *this, i ) ; }
111  inline Column column( int i ){ return Column(*this, i ) ; }
112  inline ConstColumn column( int i ) const{ return ConstColumn( *this, i ) ; }
113 
114  inline const_iterator begin() const{ return VECTOR::begin() ; }
115  inline const_iterator end() const{ return VECTOR::end() ; }
116  inline const_iterator cbegin() const{ return VECTOR::begin() ; }
117  inline const_iterator cend() const{ return VECTOR::end() ; }
118  inline iterator begin() { return VECTOR::begin() ; }
119  inline iterator end() { return VECTOR::end() ; }
120 
121  template <typename U>
122  void fill_diag( const U& u) {
124  }
125 
126  template <typename U> static Matrix diag( int size, const U& diag_value ) {
127  Matrix res(size,size) ;
128  res.fill_diag( diag_value ) ;
129  return res ;
130  }
131 
132  inline Proxy operator[]( R_xlen_t i ) {
133  return static_cast< Vector<RTYPE>* >( this )->operator[]( i ) ;
134  }
135  inline const_Proxy operator[]( R_xlen_t i ) const {
136  return static_cast< const Vector<RTYPE>* >( this )->operator[]( i ) ;
137  }
138 
139  inline Proxy operator()( const size_t& i, const size_t& j) {
140  return static_cast< Vector<RTYPE>* >( this )->operator[]( offset( i, j ) ) ;
141  }
142  inline const_Proxy operator()( const size_t& i, const size_t& j) const {
143  return static_cast< const Vector<RTYPE>* >( this )->operator[]( offset( i, j ) ) ;
144  }
145 
146  inline Proxy at( const size_t& i, const size_t& j) {
147  return static_cast< Vector<RTYPE>* >( this )->operator()( i, j ) ;
148  }
149  inline const_Proxy at( const size_t& i, const size_t& j) const {
150  return static_cast< const Vector<RTYPE>* >( this )->operator()( i, j ) ;
151  }
152 
153  inline Row operator()( int i, internal::NamedPlaceHolder ) {
154  return Row( *this, i ) ;
155  }
156  inline ConstRow operator()( int i, internal::NamedPlaceHolder ) const {
157  return ConstRow( *this, i ) ;
158  }
159  inline Column operator()( internal::NamedPlaceHolder, int i ) {
160  return Column( *this, i ) ;
161  }
162  inline ConstColumn operator()( internal::NamedPlaceHolder, int i ) const {
163  return ConstColumn( *this, i ) ;
164  }
165  inline Sub operator()( const Range& row_range, const Range& col_range) {
166  return Sub( const_cast<Matrix&>(*this), row_range, col_range ) ;
167  }
168  inline Sub operator()( internal::NamedPlaceHolder, const Range& col_range) {
169  return Sub( const_cast<Matrix&>(*this), Range(0,nrow()-1) , col_range ) ;
170  }
171  inline Sub operator()( const Range& row_range, internal::NamedPlaceHolder ) {
172  return Sub( const_cast<Matrix&>(*this), row_range, Range(0,ncol()-1) ) ;
173  }
174 
175 private:
176  inline R_xlen_t offset(const int i, const int j) const { return i + static_cast<R_xlen_t>(nrows) * j ; }
177 
178  template <typename U>
180  Shield<SEXP> elem( converter_type::get( u ) );
181 
182  R_xlen_t bounds = std::min(Matrix::nrow(), Matrix::ncol());
183  for (R_xlen_t i = 0; i < bounds; ++i) {
184  (*this)(i, i) = elem;
185  }
186  }
187 
188  template <typename U>
190  stored_type elem = converter_type::get( u );
191 
192  R_xlen_t bounds = std::min(Matrix::nrow(), Matrix::ncol());
193 
194  for (R_xlen_t i = 0; i < bounds; ++i) {
195  (*this)(i, i) = elem;
196  }
197  }
198 
199  template <bool NA, typename MAT>
200  void import_matrix_expression( const MatrixBase<RTYPE,NA,MAT>& other, int nr, int nc ) {
201  iterator start = VECTOR::begin() ;
202  for( int j=0; j<nc; j++){
203  for( int i=0; i<nr; i++, ++start){
204  *start = other(i,j) ;
205  }
206  }
207  }
208 
209 };
210 
212  return internal::DimNameProxy(x, 0);
213 }
214 
216  return internal::DimNameProxy(x, 1);
217 }
218 
219 template<template <class> class StoragePolicy >
220 inline std::ostream &operator<<(std::ostream & s, const Matrix<REALSXP, StoragePolicy> & rhs) {
221  typedef Matrix<REALSXP, StoragePolicy> MATRIX;
222 
223  std::ios::fmtflags flags = s.flags();
224  s.unsetf(std::ios::floatfield);
225  std::streamsize precision = s.precision();
226 
227  const int rows = rhs.rows();
228 
229  for (int i = 0; i < rows; ++i) {
230  typename MATRIX::Row row = const_cast<MATRIX &>(rhs).row(i);
231 
232  typename MATRIX::Row::iterator j = row.begin();
233  typename MATRIX::Row::iterator jend = row.end();
234 
235  if (j != jend) {
236  s << std::showpoint << std::setw(precision + 1) << (*j);
237  j++;
238 
239  for ( ; j != jend; j++) {
240  s << " " << std::showpoint << std::setw(precision + 1) << (*j);
241  }
242  }
243 
244  s << std::endl;
245  }
246 
247  s.flags(flags);
248  return s;
249 }
250 
251 #ifndef RCPP_NO_SUGAR
252 #define RCPP_GENERATE_MATRIX_SCALAR_OPERATOR(__OPERATOR__) \
253  template <int RTYPE, template <class> class StoragePolicy, typename T > \
254  inline typename traits::enable_if< traits::is_convertible< typename traits::remove_const_and_reference< T >::type, \
255  typename Matrix<RTYPE, StoragePolicy>::stored_type >::value, Matrix<RTYPE, StoragePolicy> >::type \
256  operator __OPERATOR__ (const Matrix<RTYPE, StoragePolicy> &lhs, const T &rhs) { \
257  Vector<RTYPE, StoragePolicy> v = static_cast<const Vector<RTYPE, StoragePolicy> &>(lhs) __OPERATOR__ rhs; \
258  v.attr("dim") = Vector<INTSXP>::create(lhs.nrow(), lhs.ncol()); \
259  return as< Matrix<RTYPE, StoragePolicy> >(v); \
260  }
261 
266 
267 #undef RCPP_GENERATE_MATRIX_SCALAR_OPERATOR
268 
269 #define RCPP_GENERATE_SCALAR_MATRIX_OPERATOR(__OPERATOR__) \
270  template <int RTYPE, template <class> class StoragePolicy, typename T > \
271  inline typename traits::enable_if< traits::is_convertible< typename traits::remove_const_and_reference< T >::type, \
272  typename Matrix<RTYPE, StoragePolicy>::stored_type >::value, Matrix<RTYPE, StoragePolicy> >::type \
273  operator __OPERATOR__ (const T &lhs, const Matrix<RTYPE, StoragePolicy> &rhs) { \
274  Vector<RTYPE, StoragePolicy> v = static_cast<const Vector<RTYPE, StoragePolicy> &>(rhs); \
275  v = lhs __OPERATOR__ v; \
276  v.attr("dim") = Vector<INTSXP>::create(rhs.nrow(), rhs.ncol()); \
277  return as< Matrix<RTYPE, StoragePolicy> >(v); \
278  }
279 
284 
285 #undef RCPP_GENERATE_SCALAR_MATRIX_OPERATOR
286 #endif
287 
288 template<template <class> class StoragePolicy >
289 inline std::ostream &operator<<(std::ostream & s, const Matrix<INTSXP, StoragePolicy> & rhs) {
290  typedef Matrix<INTSXP, StoragePolicy> MATRIX;
292 
293  std::ios::fmtflags flags = s.flags();
294 
295  s << std::dec;
296 
299 
300  typename VECTOR::iterator j = static_cast<VECTOR &>(const_cast<MATRIX &>(rhs)).begin();
301  typename VECTOR::iterator jend = static_cast<VECTOR &>(const_cast<MATRIX &>(rhs)).end();
302 
303  for ( ; j != jend; ++j) {
304  if (*j < min) {
305  min = *j;
306  }
307 
308  if (*j > max) {
309  max = *j;
310  }
311  }
312 
313  int digitsMax = (max >= 0) ? 0 : 1;
314  int digitsMin = (min >= 0) ? 0 : 1;
315 
316  while (min != 0)
317  {
318  ++digitsMin;
319  min /= 10;
320  }
321 
322  while (max != 0)
323  {
324  ++digitsMax;
325  max /= 10;
326  }
327 
328  int digits = std::max(digitsMin, digitsMax);
329 
330  const int rows = rhs.rows();
331 
332  for (int i = 0; i < rows; ++i) {
333  typename MATRIX::Row row = const_cast<MATRIX &>(rhs).row(i);
334 
335  typename MATRIX::Row::iterator j = row.begin();
336  typename MATRIX::Row::iterator jend = row.end();
337 
338  if (j != jend) {
339  s << std::setw(digits) << (*j);
340  ++j;
341 
342  for ( ; j != jend; ++j) {
343  s << " " << std::setw(digits) << (*j);
344  }
345  }
346 
347  s << std::endl;
348  }
349 
350  s.flags(flags);
351  return s;
352 }
353 
354 template<template <class> class StoragePolicy >
355 inline std::ostream &operator<<(std::ostream & s, const Matrix<STRSXP, StoragePolicy> & rhs) {
356  typedef Matrix<STRSXP, StoragePolicy> MATRIX;
357 
358  const int rows = rhs.rows();
359 
360  for (int i = 0; i < rows; ++i) {
361  typename MATRIX::Row row = const_cast<MATRIX &>(rhs).row(i);
362 
363  typename MATRIX::Row::iterator j = row.begin();
364  typename MATRIX::Row::iterator jend = row.end();
365 
366  if (j != jend) {
367  s << "\"" << (*j) << "\"";
368  j++;
369 
370  for ( ; j != jend; j++) {
371  s << " \"" << (*j) << "\"";
372  }
373  }
374 
375  s << std::endl;
376  }
377 
378  return s;
379 }
380 
381 template<int RTYPE, template <class> class StoragePolicy >
382 inline std::ostream &operator<<(std::ostream & s, const Matrix<RTYPE, StoragePolicy> & rhs) {
383  typedef Matrix<RTYPE, StoragePolicy> MATRIX;
384 
385  const int rows = rhs.rows();
386 
387  for (int i = 0; i < rows; ++i) {
388  typename MATRIX::Row row = const_cast<MATRIX &>(rhs).row(i);
389 
390  typename MATRIX::Row::iterator j = row.begin();
391  typename MATRIX::Row::iterator jend = row.end();
392 
393  if (j != jend) {
394  s << (*j);
395  j++;
396 
397  for ( ; j != jend; j++) {
398  s << (*j);
399  }
400  }
401 
402  s << std::endl;
403  }
404 
405  return s;
406 }
407 
408 template<int RTYPE, template <class> class StoragePolicy >
410  typedef Matrix<RTYPE, StoragePolicy> MATRIX;
412 
413  Vector<INTSXP, StoragePolicy> dims = ::Rf_getAttrib(x, R_DimSymbol);
414  int nrow = dims[0], ncol = dims[1];
415  MATRIX r(Dimension(ncol, nrow)); // new Matrix with reversed dimension
416  R_xlen_t len = XLENGTH(x), len2 = XLENGTH(x)-1;
417 
418  // similar approach as in R: fill by in column, "accessing row-wise"
419  VECTOR s = VECTOR(r.get__());
420  for (R_xlen_t i = 0, j = 0; i < len; i++, j += nrow) {
421  if (j > len2) j -= len2;
422  s[i] = x[j];
423  }
424 
425  // there must be a simpler, more C++-ish way for this ...
426  SEXP dimNames = Rf_getAttrib(x, R_DimNamesSymbol);
427  if (!Rf_isNull(dimNames)) {
428  // do we need dimnamesnames ?
429  Shield<SEXP> newDimNames(Rf_allocVector(VECSXP, 2));
430  SET_VECTOR_ELT(newDimNames, 0, VECTOR_ELT(dimNames, 1));
431  SET_VECTOR_ELT(newDimNames, 1, VECTOR_ELT(dimNames, 0));
432  Rf_setAttrib(r, R_DimNamesSymbol, newDimNames);
433  }
434  return r;
435 }
436 
437 template<template <class> class StoragePolicy>
439  return tranpose_impl<REALSXP, StoragePolicy>(x);
440 }
441 
442 template<template <class> class StoragePolicy>
444  return tranpose_impl<INTSXP, StoragePolicy>(x);
445 }
446 
447 template<template <class> class StoragePolicy>
449  return tranpose_impl<STRSXP, StoragePolicy>(x);
450 }
451 
452 }
453 
454 #endif
Matrix< RTYPE, StoragePolicy > tranpose_impl(const Matrix< RTYPE, StoragePolicy > &x)
Definition: Matrix.h:409
traits::enable_if< helpers::decays_to_ctype< typename std::iterator_traits< InputIterator >::value_type >::value, typename helpers::ctype< typename std::iterator_traits< InputIterator >::value_type >::type >::type min(InputIterator begin, InputIterator end)
Definition: algorithm.h:370
int ncol() const
Definition: no_init.h:55
ConstColumn column(int i) const
Definition: Matrix.h:112
Matrix(const Dimension &dims)
Definition: Matrix.h:55
void init()
Definition: Vector.h:615
traits::r_vector_const_iterator< RTYPE >::type const_iterator
Definition: Vector.h:48
const_iterator cbegin() const
Definition: Matrix.h:116
VECTOR::const_iterator const_iterator
Definition: Matrix.h:45
VECTOR::Proxy Proxy
Definition: Matrix.h:48
int * dims() const
Definition: Vector.h:611
#define RCPP_GENERATE_SCALAR_MATRIX_OPERATOR(__OPERATOR__)
Definition: Matrix.h:269
#define RCPP_GENERATE_MATRIX_SCALAR_OPERATOR(__OPERATOR__)
Definition: Matrix.h:252
SubMatrix< RTYPE > Sub
Definition: Matrix.h:40
sugar::Min< RTYPE, NA, T > min(const VectorBase< RTYPE, NA, T > &x)
Definition: min.h:80
traits::r_vector_iterator< RTYPE >::type iterator
Definition: Vector.h:47
ConstRow row(int i) const
Definition: Matrix.h:110
Sub operator()(const Range &row_range, internal::NamedPlaceHolder)
Definition: Matrix.h:171
Sub operator()(internal::NamedPlaceHolder, const Range &col_range)
Definition: Matrix.h:168
ConstColumn operator()(internal::NamedPlaceHolder, int i) const
Definition: Matrix.h:162
internal::DimNameProxy rownames(SEXP x)
Definition: Matrix.h:211
R_xlen_t size() const
Definition: Vector.h:274
int nrow() const
Definition: Matrix.h:99
Matrix(const Matrix &other)
Definition: Matrix.h:74
R_xlen_t offset(const int i, const int j) const
Definition: Matrix.h:176
const_iterator begin() const
Definition: Matrix.h:114
int nrows
Definition: Matrix.h:29
void fill_diag__dispatch(traits::false_type, const U &u)
Definition: Matrix.h:179
VECTOR::stored_type stored_type
Definition: Matrix.h:47
Matrix< REALSXP, StoragePolicy > transpose(const Matrix< REALSXP, StoragePolicy > &x)
Definition: Matrix.h:438
int rows() const
Definition: Matrix.h:105
Sub operator()(const Range &row_range, const Range &col_range)
Definition: Matrix.h:165
void import_matrix_expression(const MatrixBase< RTYPE, NA, MAT > &other, int nr, int nc)
Definition: Matrix.h:200
iterator begin()
Definition: MatrixRow.h:175
MatrixColumn< RTYPE > Column
Definition: Matrix.h:38
traits::storage_type< RTYPE >::type stored_type
Definition: Vector.h:51
Column column(int i)
Definition: Matrix.h:111
SEXP r_cast(SEXP x)
Definition: r_cast.h:151
Matrix(const MatrixBase< RTYPE, NA, MAT > &other)
Definition: Matrix.h:77
iterator begin()
Definition: Matrix.h:118
const_iterator cend() const
Definition: Matrix.h:117
const_iterator end() const
Definition: Matrix.h:115
Matrix(const int &nrows_, const int &ncols, Iterator start)
Definition: Matrix.h:64
ConstMatrixColumn< RTYPE > ConstColumn
Definition: Matrix.h:39
void fill_diag(const U &u)
Definition: Matrix.h:122
sugar::Max< RTYPE, NA, T > max(const VectorBase< RTYPE, NA, T > &x)
Definition: max.h:81
traits::enable_if< helpers::decays_to_ctype< typename std::iterator_traits< InputIterator >::value_type >::value, typename helpers::ctype< typename std::iterator_traits< InputIterator >::value_type >::type >::type max(InputIterator begin, InputIterator end)
Definition: algorithm.h:324
const_Proxy operator[](R_xlen_t i) const
Definition: Matrix.h:135
ConstMatrixRow< RTYPE > ConstRow
Definition: Matrix.h:37
iterator end()
Definition: Vector.h:333
Matrix(const no_init_matrix &obj)
Definition: Matrix.h:92
int ncol() const
Definition: Matrix.h:96
Row row(int i)
Definition: Matrix.h:109
traits::r_vector_proxy< RTYPE >::type Proxy
Definition: Vector.h:43
Matrix(SEXP x)
Definition: Matrix.h:53
int nrow() const
Definition: no_init.h:51
VECTOR::iterator iterator
Definition: Matrix.h:44
ConstRow operator()(int i, internal::NamedPlaceHolder) const
Definition: Matrix.h:156
Proxy at(const size_t &i, const size_t &j)
Definition: Matrix.h:146
Matrix & operator=(const Matrix &other)
Definition: Matrix.h:83
const_Proxy operator()(const size_t &i, const size_t &j) const
Definition: Matrix.h:142
MatrixRow< RTYPE > Row
Definition: Matrix.h:36
VECTOR::const_Proxy const_Proxy
Definition: Matrix.h:49
Matrix(const int &n)
Definition: Matrix.h:71
VECTOR::converter_type converter_type
Definition: Matrix.h:46
Rcpp API.
Definition: algo.h:28
iterator begin()
Definition: Vector.h:332
static Matrix diag(int size, const U &diag_value)
Definition: Matrix.h:126
const_Proxy at(const size_t &i, const size_t &j) const
Definition: Matrix.h:149
Proxy operator[](R_xlen_t i)
Definition: Matrix.h:132
int size() const
Definition: Dimension.h:56
StoragePolicy< Matrix > Storage
Definition: Matrix.h:42
Vector< RTYPE, StoragePolicy > VECTOR
Definition: Matrix.h:43
int cols() const
Definition: Matrix.h:102
internal::DimNameProxy colnames(SEXP x)
Definition: Matrix.h:215
Column operator()(internal::NamedPlaceHolder, int i)
Definition: Matrix.h:159
iterator end()
Definition: Matrix.h:119
Proxy operator()(const size_t &i, const size_t &j)
Definition: Matrix.h:139
Row operator()(int i, internal::NamedPlaceHolder)
Definition: Matrix.h:153
traits::r_vector_const_proxy< RTYPE >::type const_Proxy
Definition: Vector.h:44
static target get(const T &input)
Definition: converter.h:33
void fill_diag__dispatch(traits::true_type, const U &u)
Definition: Matrix.h:189
Matrix(const int &nrows_, const int &ncols)
Definition: Matrix.h:59