Rcpp Version 0.12.12
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 iterator begin() { return VECTOR::begin() ; }
117  inline iterator end() { return VECTOR::end() ; }
118 
119  template <typename U>
120  void fill_diag( const U& u) {
122  }
123 
124  template <typename U> static Matrix diag( int size, const U& diag_value ) {
125  Matrix res(size,size) ;
126  res.fill_diag( diag_value ) ;
127  return res ;
128  }
129 
130  inline Proxy operator[]( R_xlen_t i ) {
131  return static_cast< Vector<RTYPE>* >( this )->operator[]( i ) ;
132  }
133  inline const_Proxy operator[]( R_xlen_t i ) const {
134  return static_cast< const Vector<RTYPE>* >( this )->operator[]( i ) ;
135  }
136 
137  inline Proxy operator()( const size_t& i, const size_t& j) {
138  return static_cast< Vector<RTYPE>* >( this )->operator[]( offset( i, j ) ) ;
139  }
140  inline const_Proxy operator()( const size_t& i, const size_t& j) const {
141  return static_cast< const Vector<RTYPE>* >( this )->operator[]( offset( i, j ) ) ;
142  }
143 
144  inline Proxy at( const size_t& i, const size_t& j) {
145  return static_cast< Vector<RTYPE>* >( this )->operator()( i, j ) ;
146  }
147  inline const_Proxy at( const size_t& i, const size_t& j) const {
148  return static_cast< const Vector<RTYPE>* >( this )->operator()( i, j ) ;
149  }
150 
151  inline Row operator()( int i, internal::NamedPlaceHolder ) {
152  return Row( *this, i ) ;
153  }
154  inline ConstRow operator()( int i, internal::NamedPlaceHolder ) const {
155  return ConstRow( *this, i ) ;
156  }
157  inline Column operator()( internal::NamedPlaceHolder, int i ) {
158  return Column( *this, i ) ;
159  }
160  inline ConstColumn operator()( internal::NamedPlaceHolder, int i ) const {
161  return ConstColumn( *this, i ) ;
162  }
163  inline Sub operator()( const Range& row_range, const Range& col_range) {
164  return Sub( const_cast<Matrix&>(*this), row_range, col_range ) ;
165  }
166  inline Sub operator()( internal::NamedPlaceHolder, const Range& col_range) {
167  return Sub( const_cast<Matrix&>(*this), Range(0,nrow()-1) , col_range ) ;
168  }
169  inline Sub operator()( const Range& row_range, internal::NamedPlaceHolder ) {
170  return Sub( const_cast<Matrix&>(*this), row_range, Range(0,ncol()-1) ) ;
171  }
172 
173 private:
174  inline R_xlen_t offset(const int i, const int j) const { return i + static_cast<R_xlen_t>(nrows) * j ; }
175 
176  template <typename U>
178  Shield<SEXP> elem( converter_type::get( u ) );
179 
180  R_xlen_t bounds = std::min(Matrix::nrow(), Matrix::ncol());
181  for (R_xlen_t i = 0; i < bounds; ++i) {
182  (*this)(i, i) = elem;
183  }
184  }
185 
186  template <typename U>
188  stored_type elem = converter_type::get( u );
189 
190  R_xlen_t bounds = std::min(Matrix::nrow(), Matrix::ncol());
191 
192  for (R_xlen_t i = 0; i < bounds; ++i) {
193  (*this)(i, i) = elem;
194  }
195  }
196 
197  template <bool NA, typename MAT>
198  void import_matrix_expression( const MatrixBase<RTYPE,NA,MAT>& other, int nr, int nc ) {
199  iterator start = VECTOR::begin() ;
200  for( int j=0; j<nc; j++){
201  for( int i=0; i<nr; i++, ++start){
202  *start = other(i,j) ;
203  }
204  }
205  }
206 
207 };
208 
210  return internal::DimNameProxy(x, 0);
211 }
212 
214  return internal::DimNameProxy(x, 1);
215 }
216 
217 template<template <class> class StoragePolicy >
218 inline std::ostream &operator<<(std::ostream & s, const Matrix<REALSXP, StoragePolicy> & rhs) {
219  typedef Matrix<REALSXP, StoragePolicy> MATRIX;
220 
221  std::ios::fmtflags flags = s.flags();
222  s.unsetf(std::ios::floatfield);
223  std::streamsize precision = s.precision();
224 
225  const int rows = rhs.rows();
226 
227  for (int i = 0; i < rows; ++i) {
228  typename MATRIX::Row row = const_cast<MATRIX &>(rhs).row(i);
229 
230  typename MATRIX::Row::iterator j = row.begin();
231  typename MATRIX::Row::iterator jend = row.end();
232 
233  if (j != jend) {
234  s << std::showpoint << std::setw(precision + 1) << (*j);
235  j++;
236 
237  for ( ; j != jend; j++) {
238  s << " " << std::showpoint << std::setw(precision + 1) << (*j);
239  }
240  }
241 
242  s << std::endl;
243  }
244 
245  s.flags(flags);
246  return s;
247 }
248 
249 #ifndef RCPP_NO_SUGAR
250 #define RCPP_GENERATE_MATRIX_SCALAR_OPERATOR(__OPERATOR__) \
251  template <int RTYPE, template <class> class StoragePolicy, typename T > \
252  inline typename traits::enable_if< traits::is_convertible< typename traits::remove_const_and_reference< T >::type, \
253  typename Matrix<RTYPE, StoragePolicy>::stored_type >::value, Matrix<RTYPE, StoragePolicy> >::type \
254  operator __OPERATOR__ (const Matrix<RTYPE, StoragePolicy> &lhs, const T &rhs) { \
255  Vector<RTYPE, StoragePolicy> v = static_cast<const Vector<RTYPE, StoragePolicy> &>(lhs) __OPERATOR__ rhs; \
256  v.attr("dim") = Vector<INTSXP>::create(lhs.nrow(), lhs.ncol()); \
257  return as< Matrix<RTYPE, StoragePolicy> >(v); \
258  }
259 
264 
265 #undef RCPP_GENERATE_MATRIX_SCALAR_OPERATOR
266 
267 #define RCPP_GENERATE_SCALAR_MATRIX_OPERATOR(__OPERATOR__) \
268  template <int RTYPE, template <class> class StoragePolicy, typename T > \
269  inline typename traits::enable_if< traits::is_convertible< typename traits::remove_const_and_reference< T >::type, \
270  typename Matrix<RTYPE, StoragePolicy>::stored_type >::value, Matrix<RTYPE, StoragePolicy> >::type \
271  operator __OPERATOR__ (const T &lhs, const Matrix<RTYPE, StoragePolicy> &rhs) { \
272  Vector<RTYPE, StoragePolicy> v = static_cast<const Vector<RTYPE, StoragePolicy> &>(rhs); \
273  v = lhs __OPERATOR__ v; \
274  v.attr("dim") = Vector<INTSXP>::create(rhs.nrow(), rhs.ncol()); \
275  return as< Matrix<RTYPE, StoragePolicy> >(v); \
276  }
277 
282 
283 #undef RCPP_GENERATE_SCALAR_MATRIX_OPERATOR
284 #endif
285 
286 template<template <class> class StoragePolicy >
287 inline std::ostream &operator<<(std::ostream & s, const Matrix<INTSXP, StoragePolicy> & rhs) {
288  typedef Matrix<INTSXP, StoragePolicy> MATRIX;
290 
291  std::ios::fmtflags flags = s.flags();
292 
293  s << std::dec;
294 
297 
298  typename VECTOR::iterator j = static_cast<VECTOR &>(const_cast<MATRIX &>(rhs)).begin();
299  typename VECTOR::iterator jend = static_cast<VECTOR &>(const_cast<MATRIX &>(rhs)).end();
300 
301  for ( ; j != jend; ++j) {
302  if (*j < min) {
303  min = *j;
304  }
305 
306  if (*j > max) {
307  max = *j;
308  }
309  }
310 
311  int digitsMax = (max >= 0) ? 0 : 1;
312  int digitsMin = (min >= 0) ? 0 : 1;
313 
314  while (min != 0)
315  {
316  ++digitsMin;
317  min /= 10;
318  }
319 
320  while (max != 0)
321  {
322  ++digitsMax;
323  max /= 10;
324  }
325 
326  int digits = std::max(digitsMin, digitsMax);
327 
328  const int rows = rhs.rows();
329 
330  for (int i = 0; i < rows; ++i) {
331  typename MATRIX::Row row = const_cast<MATRIX &>(rhs).row(i);
332 
333  typename MATRIX::Row::iterator j = row.begin();
334  typename MATRIX::Row::iterator jend = row.end();
335 
336  if (j != jend) {
337  s << std::setw(digits) << (*j);
338  ++j;
339 
340  for ( ; j != jend; ++j) {
341  s << " " << std::setw(digits) << (*j);
342  }
343  }
344 
345  s << std::endl;
346  }
347 
348  s.flags(flags);
349  return s;
350 }
351 
352 template<template <class> class StoragePolicy >
353 inline std::ostream &operator<<(std::ostream & s, const Matrix<STRSXP, StoragePolicy> & rhs) {
354  typedef Matrix<STRSXP, StoragePolicy> MATRIX;
355 
356  const int rows = rhs.rows();
357 
358  for (int i = 0; i < rows; ++i) {
359  typename MATRIX::Row row = const_cast<MATRIX &>(rhs).row(i);
360 
361  typename MATRIX::Row::iterator j = row.begin();
362  typename MATRIX::Row::iterator jend = row.end();
363 
364  if (j != jend) {
365  s << "\"" << (*j) << "\"";
366  j++;
367 
368  for ( ; j != jend; j++) {
369  s << " \"" << (*j) << "\"";
370  }
371  }
372 
373  s << std::endl;
374  }
375 
376  return s;
377 }
378 
379 template<int RTYPE, template <class> class StoragePolicy >
380 inline std::ostream &operator<<(std::ostream & s, const Matrix<RTYPE, StoragePolicy> & rhs) {
381  typedef Matrix<RTYPE, StoragePolicy> MATRIX;
382 
383  const int rows = rhs.rows();
384 
385  for (int i = 0; i < rows; ++i) {
386  typename MATRIX::Row row = const_cast<MATRIX &>(rhs).row(i);
387 
388  typename MATRIX::Row::iterator j = row.begin();
389  typename MATRIX::Row::iterator jend = row.end();
390 
391  if (j != jend) {
392  s << (*j);
393  j++;
394 
395  for ( ; j != jend; j++) {
396  s << (*j);
397  }
398  }
399 
400  s << std::endl;
401  }
402 
403  return s;
404 }
405 
406 template<int RTYPE, template <class> class StoragePolicy >
408  typedef Matrix<RTYPE, StoragePolicy> MATRIX;
410 
411  Vector<INTSXP, StoragePolicy> dims = ::Rf_getAttrib(x, R_DimSymbol);
412  int nrow = dims[0], ncol = dims[1];
413  MATRIX r(Dimension(ncol, nrow)); // new Matrix with reversed dimension
414  R_xlen_t len = XLENGTH(x), len2 = XLENGTH(x)-1;
415 
416  // similar approach as in R: fill by in column, "accessing row-wise"
417  VECTOR s = VECTOR(r.get__());
418  for (R_xlen_t i = 0, j = 0; i < len; i++, j += nrow) {
419  if (j > len2) j -= len2;
420  s[i] = x[j];
421  }
422 
423  // there must be a simpler, more C++-ish way for this ...
424  SEXP dimNames = Rf_getAttrib(x, R_DimNamesSymbol);
425  if (!Rf_isNull(dimNames)) {
426  // do we need dimnamesnames ?
427  Shield<SEXP> newDimNames(Rf_allocVector(VECSXP, 2));
428  SET_VECTOR_ELT(newDimNames, 0, VECTOR_ELT(dimNames, 1));
429  SET_VECTOR_ELT(newDimNames, 1, VECTOR_ELT(dimNames, 0));
430  Rf_setAttrib(r, R_DimNamesSymbol, newDimNames);
431  }
432  return r;
433 }
434 
435 template<template <class> class StoragePolicy>
437  return tranpose_impl<REALSXP, StoragePolicy>(x);
438 }
439 
440 template<template <class> class StoragePolicy>
442  return tranpose_impl<INTSXP, StoragePolicy>(x);
443 }
444 
445 template<template <class> class StoragePolicy>
447  return tranpose_impl<STRSXP, StoragePolicy>(x);
448 }
449 
450 }
451 
452 #endif
Matrix< RTYPE, StoragePolicy > tranpose_impl(const Matrix< RTYPE, StoragePolicy > &x)
Definition: Matrix.h:407
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
Matrix(const Dimension &dims)
Definition: Matrix.h:55
void init()
Definition: Vector.h:613
const_iterator begin() const
Definition: Matrix.h:114
traits::r_vector_const_iterator< RTYPE >::type const_iterator
Definition: Vector.h:48
int nrow() const
Definition: Matrix.h:99
VECTOR::const_iterator const_iterator
Definition: Matrix.h:45
int * dims() const
Definition: Vector.h:609
ConstColumn operator()(internal::NamedPlaceHolder, int i) const
Definition: Matrix.h:160
VECTOR::Proxy Proxy
Definition: Matrix.h:48
#define RCPP_GENERATE_SCALAR_MATRIX_OPERATOR(__OPERATOR__)
Definition: Matrix.h:267
#define RCPP_GENERATE_MATRIX_SCALAR_OPERATOR(__OPERATOR__)
Definition: Matrix.h:250
int nrow() const
Definition: no_init.h:51
int rows() const
Definition: Matrix.h:105
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
Sub operator()(const Range &row_range, internal::NamedPlaceHolder)
Definition: Matrix.h:169
Sub operator()(internal::NamedPlaceHolder, const Range &col_range)
Definition: Matrix.h:166
const_Proxy at(const size_t &i, const size_t &j) const
Definition: Matrix.h:147
ConstRow operator()(int i, internal::NamedPlaceHolder) const
Definition: Matrix.h:154
internal::DimNameProxy rownames(SEXP x)
Definition: Matrix.h:209
Matrix(const Matrix &other)
Definition: Matrix.h:74
int cols() const
Definition: Matrix.h:102
ConstColumn column(int i) const
Definition: Matrix.h:112
int nrows
Definition: Matrix.h:29
void fill_diag__dispatch(traits::false_type, const U &u)
Definition: Matrix.h:177
VECTOR::stored_type stored_type
Definition: Matrix.h:47
Matrix< REALSXP, StoragePolicy > transpose(const Matrix< REALSXP, StoragePolicy > &x)
Definition: Matrix.h:436
Sub operator()(const Range &row_range, const Range &col_range)
Definition: Matrix.h:163
int size() const
Definition: Dimension.h:56
void import_matrix_expression(const MatrixBase< RTYPE, NA, MAT > &other, int nr, int nc)
Definition: Matrix.h:198
iterator begin()
Definition: MatrixRow.h:146
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
R_xlen_t size() const
Definition: Vector.h:274
iterator begin()
Definition: Matrix.h:116
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:120
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
ConstRow row(int i) const
Definition: Matrix.h:110
ConstMatrixRow< RTYPE > ConstRow
Definition: Matrix.h:37
iterator end()
Definition: Vector.h:333
Matrix(const no_init_matrix &obj)
Definition: Matrix.h:92
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
const_Proxy operator[](R_xlen_t i) const
Definition: Matrix.h:133
VECTOR::iterator iterator
Definition: Matrix.h:44
Proxy at(const size_t &i, const size_t &j)
Definition: Matrix.h:144
Matrix & operator=(const Matrix &other)
Definition: Matrix.h:83
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:124
Proxy operator[](R_xlen_t i)
Definition: Matrix.h:130
StoragePolicy< Matrix > Storage
Definition: Matrix.h:42
Vector< RTYPE, StoragePolicy > VECTOR
Definition: Matrix.h:43
internal::DimNameProxy colnames(SEXP x)
Definition: Matrix.h:213
Column operator()(internal::NamedPlaceHolder, int i)
Definition: Matrix.h:157
iterator end()
Definition: Matrix.h:117
Proxy operator()(const size_t &i, const size_t &j)
Definition: Matrix.h:137
Row operator()(int i, internal::NamedPlaceHolder)
Definition: Matrix.h:151
int ncol() const
Definition: no_init.h:55
traits::r_vector_const_proxy< RTYPE >::type const_Proxy
Definition: Vector.h:44
R_xlen_t offset(const int i, const int j) const
Definition: Matrix.h:174
const_iterator end() const
Definition: Matrix.h:115
static target get(const T &input)
Definition: converter.h:33
void fill_diag__dispatch(traits::true_type, const U &u)
Definition: Matrix.h:187
const_Proxy operator()(const size_t &i, const size_t &j) const
Definition: Matrix.h:140
Matrix(const int &nrows_, const int &ncols)
Definition: Matrix.h:59
int ncol() const
Definition: Matrix.h:96