1 /*
2  * Scythe Statistical Library Copyright (C) 2000-2002 Andrew D. Martin
3  * and Kevin M. Quinn; 2002-present Andrew D. Martin, Kevin M. Quinn,
4  * and Daniel Pemstein.  All Rights Reserved.
5  *
6  * This program is free software; you can redistribute it and/or
7  * modify under the terms of the GNU General Public License as
8  * published by Free Software Foundation; either version 2 of the
9  * License, or (at your option) any later version.  See the text files
10  * COPYING and LICENSE, distributed with this source code, for further
11  * information.
12  * --------------------------------------------------------------------
13  *  scythestat/matrix_random_access_iterator.h
14  *
15  * Random access iterators for the matrix class.
16  *
17  */
18 
19 /*! \file matrix_random_access_iterator.h
20  * \brief Definitions of STL-compliant random access iterators for
21  * the Matrix class.
22  *
23  * Contains definitions of const_matrix_random_access_iterator,
24  * matrix_random_access_iterator, and related operators.  See a
25  * Standard Template Library reference, such as Josuttis (1999), for a
26  * full description of the capabilities of random access iterators.
27  *
28  * These iterators are templated on the type, order and style of the
29  * Matrix they iterate over and their own order, which need not match
30  * the iterated-over matrix.  Same-order iteration over concrete
31  * matrices is extremely fast.  Cross-grain concrete and/or view
32  * iteration is slower.
33  */
34 
35 #ifndef SCYTHE_MATRIX_RANDOM_ACCESS_ITERATOR_H
36 #define SCYTHE_MATRIX_RANDOM_ACCESS_ITERATOR_H
37 
38 #include <iterator>
39 
40 #ifdef SCYTHE_COMPILE_DIRECT
41 #include "defs.h"
42 #include "error.h"
43 #include "matrix.h"
44 #else
45 #include "scythestat/defs.h"
46 #include "scythestat/error.h"
47 #include "scythestat/matrix.h"
48 #endif
49 
50 /* The const_matrix_iterator and matrix_iterator classes are
51  * essentially identical, except for the return types of the *, ->,
52  * and [] operators.  matrix_iterator extends const_matrix_iterator,
53  * overriding most of its members. */
54 
55 /* TODO Current setup uses template argument based branches to
56  * handle views and cross-grained orderings differently than simple
57  * in-order concrete matrices.  The work for this gets done at
58  * compile time, but we end with a few unused instance variables in
59  * the concrete case.  It might be better to specialize the entire
60  * class, although this will lead to a lot of code duplication.  We
61  * should bench the difference down the road and see if it is worth
62  * the maintenance hassle.
63  *
64  * At the moment this is looking like it won't be worth it.
65  * Iterator-based operations on concretes provide comparable
66  * performance to element-access based routines in previous versions
67  * of the library, indicating little performance penalty.
68  */
69 
70 namespace scythe {
71 	/* convenience typedefs */
72   namespace { // local to this file
73     typedef unsigned int uint;
74   }
75 
76 	/* forward declaration of the matrix class */
77 	template <typename T_type, matrix_order ORDER, matrix_style STYLE>
78 	class Matrix;
79 
80   /*! \brief An STL-compliant const random access iterator for Matrix.
81    *
82    * Provides random access iteration over const Matrix objects.  See
83    * Josuttis (1999), or some other STL reference, for a full
84    * description of the random access iterator interface.
85    *
86    * \see Matrix
87    * \see matrix_random_access_iterator
88    * \see const_matrix_forward_iterator
89    * \see matrix_forward_iterator
90    * \see const_matrix_bidirectional_iterator
91    * \see matrix_bidirectional_iterator
92    */
93 
94   template <typename T_type, matrix_order ORDER, matrix_order M_ORDER,
95             matrix_style M_STYLE>
96   class const_matrix_random_access_iterator
97     : public std::iterator<std::random_access_iterator_tag, T_type>
98   {
99 		public:
100 			/**** TYPEDEFS ***/
101 			typedef const_matrix_random_access_iterator<T_type, ORDER,
102               M_ORDER, M_STYLE> self;
103 
104 			/* These are a little formal, but useful */
105 			typedef typename std::iterator_traits<self>::value_type
106 				value_type;
107 			typedef typename std::iterator_traits<self>::iterator_category
108 				iterator_category;
109 			typedef typename std::iterator_traits<self>::difference_type
110 				difference_type;
111 			typedef typename std::iterator_traits<self>::pointer pointer;
112 			typedef typename std::iterator_traits<self>::reference reference;
113 
114 
115 			/**** CONSTRUCTORS ****/
116 
117 			/* Default constructor */
const_matrix_random_access_iterator()118 			const_matrix_random_access_iterator ()
119 			{}
120 
121 			/* Standard constructor */
const_matrix_random_access_iterator(const Matrix<value_type,M_ORDER,M_STYLE> & M)122 			const_matrix_random_access_iterator
123         ( const Matrix<value_type, M_ORDER, M_STYLE> &M)
124         : start_ (M.getArray())
125       {
126         SCYTHE_CHECK_30 (start_ == 0, scythe_null_error,
127             "Requesting iterator to NULL matrix");
128         pos_ = start_;
129 
130         /* The basic story is: when M_STYLE == Concrete and ORDER ==
131          * M_ORDER, we only need pos_ and start_ and iteration will be
132          * as fast as possible.  All other types of iteration need
133          * more variables to keep track of things and are slower.
134          */
135 
136         if (M_STYLE != Concrete || M_ORDER != ORDER) {
137           offset_ = 0;
138 
139           if (ORDER == Col) {
140             lead_length_ = M.rows();
141             lead_inc_ = M.rowstride();
142             trail_inc_ = M.colstride();
143           } else {
144             lead_length_ = M.cols();
145             lead_inc_ = M.colstride();
146             trail_inc_ = M.rowstride();
147           }
148           jump_ = trail_inc_ + (1 - lead_length_) * lead_inc_;
149         }
150 
151 #if SCYTHE_DEBUG > 2
152 				size_ = M.size();
153 #endif
154       }
155 
156       /* Copy constructor */
const_matrix_random_access_iterator(const self & mi)157       const_matrix_random_access_iterator (const self &mi)
158         : start_ (mi.start_),
159           pos_ (mi.pos_)
160       {
161         if (M_STYLE != Concrete || M_ORDER != ORDER) {
162           offset_ = mi.offset_;
163           lead_length_ = mi.lead_length_;
164           lead_inc_ = mi.lead_inc_;
165           trail_inc_ = mi.trail_inc_;
166           jump_ = mi.jump_;
167         }
168 #if SCYTHE_DEBUG > 2
169 				size_ = mi.size_;
170 #endif
171       }
172 
173       /**** FORWARD ITERATOR FACILITIES ****/
174 
175       inline self& operator= (const self& mi)
176       {
177         start_ = mi.start_;
178         pos_ = mi.pos_;
179 
180         if (M_STYLE != Concrete || M_ORDER != ORDER) {
181           offset_ = mi.offset_;
182           lead_length_ = mi.lead_length_;
183           lead_inc_ = mi.lead_inc_;
184           trail_inc_ = mi.trail_inc_;
185           jump_ = mi.jump_;
186         }
187 #if SCYTHE_DEBUG > 2
188 				size_ = mi.size_;
189 #endif
190 
191         return *this;
192       }
193 
194       inline const reference operator* () const
195       {
196 				SCYTHE_ITER_CHECK_BOUNDS();
197         return *pos_;
198       }
199 
200       inline const pointer operator-> () const
201       {
202 				SCYTHE_ITER_CHECK_BOUNDS();
203         return pos_;
204       }
205 
206       inline self& operator++ ()
207       {
208         if (M_STYLE == Concrete && ORDER == M_ORDER)
209           ++pos_;
210         else if (++offset_ % lead_length_ == 0)
211           pos_ += jump_;
212         else
213           pos_ += lead_inc_;
214 
215         return *this;
216       }
217 
218       inline self operator++ (int)
219       {
220         self tmp = *this;
221         ++(*this);
222         return tmp;
223       }
224 
225       /* == is only defined for iterators of the same template type
226        * that point to the same matrix.  Behavior for any other
227        * comparison is undefined.
228        *
229        * Note that we have to be careful about iterator comparisons
230        * when working with views and cross-grain iterators.
231        * Specifically, we always have to rely on the offset value.
232        * Obviously, with <> checks pos_ can jump all over the place in
233        * cross-grain iterators, but also end iterators point to the
234        * value after the last in the matrix.  In some cases, the
235        * equation in += and -= will actually put pos_ inside the
236        * matrix (often in an early position) in this case.
237        */
238       inline bool operator== (const self& x) const
239       {
240         if (M_STYLE == Concrete && ORDER == M_ORDER) {
241           return pos_ == x.pos_;
242         } else {
243           return offset_ == x.offset_;
244         }
245       }
246 
247       /* Again, != is only officially defined for iterators over the
248        * same matrix although the test will be trivially true for
249        * matrices that don't view the same data, by implementation.
250        */
251       inline bool operator!= (const self &x) const
252       {
253         return !(*this == x);
254       }
255 
256       /**** BIDIRECTIONAL ITERATOR FACILITIES ****/
257 
258       inline self& operator-- ()
259       {
260         if (M_STYLE == Concrete && ORDER == M_ORDER)
261           --pos_;
262         else if (offset_-- % lead_length_ == 0)
263           pos_ -= jump_;
264         else
265           pos_ -= lead_inc_;
266 
267         return *this;
268       }
269 
270       inline self operator-- (int)
271       {
272         self tmp = *this;
273         --(*this);
274         return tmp;
275       }
276 
277       /**** RANDOM ACCESS ITERATOR FACILITIES ****/
278 
279       inline const reference operator[] (difference_type n) const
280       {
281         if (M_STYLE == Concrete && ORDER == M_ORDER) {
282 					SCYTHE_ITER_CHECK_OFFSET_BOUNDS(start_ + n);
283           return *(start_ + n);
284         } else {
285           uint trailing = n / lead_length_;
286           uint leading = n % lead_length_;
287 
288           T_type* place = start_ + leading * lead_inc_
289                                  + trailing * trail_inc_;
290 
291 					SCYTHE_ITER_CHECK_POINTER_BOUNDS(place);
292           return *place;
293         }
294       }
295 
296       inline self& operator+= (difference_type n)
297       {
298         if (M_STYLE == Concrete && ORDER == M_ORDER) {
299           pos_ += n;
300         } else {
301           offset_ += n;
302           uint trailing = offset_ / lead_length_;
303           uint leading = offset_ % lead_length_;
304 
305           pos_ = start_ + leading * lead_inc_
306                         + trailing * trail_inc_;
307         }
308 
309         return *this;
310       }
311 
312       inline self& operator-= (difference_type n)
313       {
314         if (M_STYLE == Concrete && ORDER == M_ORDER) {
315           pos_ -= n;
316         } else {
317           offset_ -= n;
318           uint trailing = offset_ / lead_length_;
319           uint leading = offset_ % lead_length_;
320 
321           pos_ = start_ + leading * lead_inc_
322                         + trailing * trail_inc_;
323         }
324 
325         return *this;
326       }
327 
328       /* + and - difference operators are outside the class */
329 
330       inline difference_type operator- (const self& x) const
331       {
332         if (M_STYLE == Concrete && ORDER == M_ORDER) {
333           return pos_ - x.pos_;
334         } else {
335           return offset_ - x.offset_;
336         }
337       }
338 
339       inline difference_type operator< (const self& x) const
340       {
341         if (M_STYLE == Concrete && ORDER == M_ORDER) {
342           return pos_ < x.pos_;
343         } else {
344           return offset_ < x.offset_;
345         }
346       }
347 
348       inline difference_type operator> (const self& x) const
349       {
350         if (M_STYLE == Concrete && ORDER == M_ORDER) {
351           return pos_ > x.pos_;
352         } else {
353           return offset_ > x.offset_;
354         }
355       }
356 
357       inline difference_type operator<= (const self& x) const
358       {
359         if (M_STYLE == Concrete && ORDER == M_ORDER) {
360           return pos_ <= x.pos_;
361         } else {
362           return offset_ <= x.offset_;
363         }
364       }
365 
366       inline difference_type operator>= (const self& x) const
367       {
368         if (M_STYLE == Concrete && ORDER == M_ORDER) {
369           return pos_ >= x.pos_;
370         } else {
371           return offset_ >= x.offset_;
372         }
373       }
374 
375     protected:
376 
377       /**** INSTANCE VARIABLES ****/
378       T_type* start_; // pointer to beginning of data array
379       T_type* pos_;         // pointer to current position in array
380 
381       uint offset_;  // Logical offset into matrix
382 
383       // TODO Some of these can probably be uints
384       int lead_length_;  // Logical length of leading dimension
385       int lead_inc_;  // Memory distance between vectors in ldim
386       int trail_inc_; // Memory distance between vectors in tdim
387       int jump_; // Memory distance between end of one ldim vector and
388                  // begin of next
389 
390 			// Size variable for range checking
391 #if SCYTHE_DEBUG > 2
392 			uint size_;  // Logical matrix size
393 #endif
394 
395  };
396 
397   /*! \brief An STL-compliant random access iterator for Matrix.
398    *
399    * Provides random access iteration over Matrix objects.  See
400    * Josuttis (1999), or some other STL reference, for a full
401    * description of the random access iterator interface.
402    *
403    * \see Matrix
404    * \see const_matrix_random_access_iterator
405    * \see const_matrix_forward_iterator
406    * \see matrix_forward_iterator
407    * \see const_matrix_bidirectional_iterator
408    * \see matrix_bidirectional_iterator
409    */
410 	template <typename T_type, matrix_order ORDER, matrix_order M_ORDER,
411             matrix_style M_STYLE>
412 	class matrix_random_access_iterator
413 		: public const_matrix_random_access_iterator<T_type, ORDER,
414                                                  M_ORDER, M_STYLE>
415 	{
416 			/**** TYPEDEFS ***/
417 			typedef matrix_random_access_iterator<T_type, ORDER, M_ORDER,
418                                             M_STYLE> self;
419 			typedef const_matrix_random_access_iterator<T_type, ORDER,
420                                                  M_ORDER, M_STYLE> Base;
421 
422 		public:
423 			/* These are a little formal, but useful */
424 			typedef typename std::iterator_traits<Base>::value_type
425 				value_type;
426 			typedef typename std::iterator_traits<Base>::iterator_category
427 				iterator_category;
428 			typedef typename std::iterator_traits<Base>::difference_type
429 				difference_type;
430 			typedef typename std::iterator_traits<Base>::pointer pointer;
431 			typedef typename std::iterator_traits<Base>::reference reference;
432 
433 
434 			/**** CONSTRUCTORS ****/
435 
436 			/* Default constructor */
matrix_random_access_iterator()437 			matrix_random_access_iterator ()
438 				: Base ()
439 			{}
440 
441 			/* Standard constructor */
matrix_random_access_iterator(const Matrix<value_type,M_ORDER,M_STYLE> & M)442 			matrix_random_access_iterator (const Matrix<value_type, M_ORDER,
443                                                   M_STYLE> &M)
444 				:	Base(M)
445 			{}
446 
447       /* Copy constructor */
matrix_random_access_iterator(const self & mi)448 			matrix_random_access_iterator (const self &mi)
449 				:	Base (mi)
450 			{}
451 
452 			/**** FORWARD ITERATOR FACILITIES ****/
453 
454 			/* We have to override a lot of these to get return values
455 			 * right.*/
456       inline self& operator= (const self& mi)
457       {
458         start_ = mi.start_;
459         pos_ = mi.pos_;
460 
461         if (M_STYLE != Concrete || M_ORDER != ORDER) {
462           offset_ = mi.offset_;
463           lead_length_ = mi.lead_length_;
464           lead_inc_ = mi.lead_inc_;
465           trail_inc_ = mi.trail_inc_;
466           jump_ = mi.jump_;
467         }
468 #if SCYTHE_DEBUG > 2
469 				size_ = mi.size_;
470 #endif
471 
472         return *this;
473       }
474 
475 			inline reference operator* () const
476 			{
477 				SCYTHE_ITER_CHECK_BOUNDS();
478 				return *pos_;
479 			}
480 
481 			inline pointer operator-> () const
482 			{
483 				SCYTHE_ITER_CHECK_BOUNDS();
484 				return pos_;
485 			}
486 
487 			inline self& operator++ ()
488 			{
489 				Base::operator++();
490 				return *this;
491 			}
492 
493 			inline self operator++ (int)
494 			{
495 				self tmp = *this;
496 				++(*this);
497 				return tmp;
498 			}
499 
500 			/**** BIDIRECTIONAL ITERATOR FACILITIES ****/
501 
502 			inline self& operator-- ()
503 			{
504 				Base::operator--();
505 				return *this;
506 			}
507 
508 			inline self operator-- (int)
509 			{
510 				self tmp = *this;
511 				--(*this);
512 				return tmp;
513 			}
514 
515 			/**** RANDOM ACCESS ITERATOR FACILITIES ****/
516 
517       inline reference operator[] (difference_type n) const
518       {
519         if (M_STYLE == Concrete && ORDER == M_ORDER) {
520 					SCYTHE_ITER_CHECK_POINTER_BOUNDS(start_ + n);
521           return *(start_ + n);
522         } else {
523           uint trailing = n / lead_length_;
524           uint leading = n % lead_length_;
525 
526           T_type* place = start_ + leading * lead_inc_
527                                  + trailing * trail_inc_;
528 
529 					SCYTHE_ITER_CHECK_POINTER_BOUNDS(place);
530           return *place;
531         }
532       }
533 
534 			inline self& operator+= (difference_type n)
535 			{
536 				Base::operator+=(n);
537 				return *this;
538 			}
539 
540 			inline self& operator-= (difference_type n)
541 			{
542 				Base::operator-= (n);
543 				return *this;
544 			}
545 
546 			/* +  and - difference_type operators are outside the class */
547 
548 		private:
549 			/* Get handles to base members.  It boggles the mind */
550 			using Base::start_;
551 			using Base::pos_;
552       using Base::offset_;
553       using Base::lead_length_;
554       using Base::lead_inc_;
555       using Base::trail_inc_;
556       using Base::jump_;
557 #if SCYTHE_DEBUG > 2
558 			using Base::size_;
559 #endif
560 	};
561 
562 	template <class T_type, matrix_order ORDER, matrix_order M_ORDER,
563             matrix_style STYLE>
564 	inline
565   const_matrix_random_access_iterator<T_type, ORDER, M_ORDER, STYLE>
566 	operator+ (const_matrix_random_access_iterator<T_type, ORDER, M_ORDER,                                                 STYLE> x, int n)
567 	{
568 		x += n;
569 		return x;
570 	}
571 
572 	template <class T_type, matrix_order ORDER, matrix_order M_ORDER,
573             matrix_style STYLE>
574 	inline
575   const_matrix_random_access_iterator<T_type, ORDER, M_ORDER, STYLE>
576 	operator+ (int n, const_matrix_random_access_iterator<T_type, ORDER,
577                                                         M_ORDER,
578                                                         STYLE> x)
579 	{
580 		x += n;
581 		return x;
582 	}
583 
584 	template <class T_type, matrix_order ORDER, matrix_order M_ORDER,
585             matrix_style STYLE>
586 	inline
587   const_matrix_random_access_iterator<T_type, ORDER, M_ORDER, STYLE>
588 	operator- (const_matrix_random_access_iterator<T_type, ORDER, M_ORDER,                                                 STYLE> x, int n)
589 	{
590 		x -= n;
591 		return x;
592 	}
593 
594 	template <class T_type, matrix_order ORDER, matrix_order M_ORDER,
595             matrix_style STYLE>
596 	inline matrix_random_access_iterator<T_type, ORDER, M_ORDER, STYLE>
597 	operator+ (matrix_random_access_iterator<T_type, ORDER, M_ORDER,
598                                            STYLE> x, int n)
599 	{
600 		x += n;
601 		return x;
602 	}
603 
604 	template <class T_type, matrix_order ORDER, matrix_order M_ORDER,
605             matrix_style STYLE>
606 	inline matrix_random_access_iterator<T_type, ORDER, M_ORDER, STYLE>
607 	operator+ (int n, matrix_random_access_iterator<T_type, ORDER,
608                                                   M_ORDER, STYLE> x)
609 	{
610 		x += n;
611 		return x;
612 	}
613 
614 	template <class T_type, matrix_order ORDER, matrix_order M_ORDER,
615             matrix_style STYLE>
616 	inline matrix_random_access_iterator<T_type, ORDER, M_ORDER, STYLE>
617 	operator- (matrix_random_access_iterator<T_type, ORDER, M_ORDER,
618                                            STYLE> x, int n)
619 	{
620 		x -= n;
621 		return x;
622 	}
623 
624 } // namespace scythe
625 
626 #endif /* SCYTHE_MATRIX_ITERATOR_H */
627