reduce Eigen headers visibility
[scilab.git] / scilab / modules / ast / includes / types / matrixiterator.hxx
1 /*
2 *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 *  Copyright (C) 2008-2008 - DIGITEO - Bernard HUGUENEY
4 *
5 *  This file must be used under the terms of the CeCILL.
6 *  This source file is licensed as described in the file COPYING, which
7 *  you should have received as part of this distribution.  The terms
8 *  are also available at
9 *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 *
11 */
12 #ifndef MATRIXITERATORS_HXX
13 #define MATRIXITERATORS_HXX
14
15 #include <complex>
16 #include <utility>
17 #include <iterator>
18
19 #include <Eigen/Sparse>
20 #include "double.hxx"
21 #include "sparse.hxx"
22
23
24 /*
25   In order to reuse code for the various Matrix Classes, we need some uniform API to access elements.
26   We cannot use runtime polymorphism (with dynamic dispatching) because of the runtime cost so we have to
27   use compile-time polymorphism with templates.
28
29   The provided free function templates get<>() and set<>() provide such an uniform API.
30
31   In order to perform element-wise operations on Matrices (copy, partial or total assignments, etc.),
32   we provide an iterator. To enable reading (with get<>()) or writiting (with set<>()) we provide an Accessor<> proxy.
33
34   As is it common to iterate over a sub-matrix according to indices given by a Scilab variable (Double) we provide
35   an iterator created from such a variable (IteratorFromVar).
36  */
37 template<typename T>
38 struct UndefinedAccessorForType {};
39
40 /**
41    This free function overloads perform read access into a 2D container, using 0-based indices.
42    @param s the 2D structure used to fetch a value of type V.
43    @param r : the row (0 based)
44    @param c : the column (0 based)
45    @return : the value of type V at row r and column c of structure s
46 */
47
48 template<typename V, typename S> V get(S SPARSE_CONST&, int, int)
49 {
50     return UndefinedAccessorForType<S>();
51 }
52
53 template<> double get(types::Double SPARSE_CONST& d, int r, int c)
54 {
55     return d.getReal(r, c);
56 }
57 template<> std::complex<double> get(types::Double SPARSE_CONST& d, int r, int c)
58 {
59     return std::complex<double>(d.getReal(r, c), d.getImg(r, c));
60 }
61
62 template<> bool get(types::Bool SPARSE_CONST& d, int r, int c)
63 {
64     return d.get(r, c) == 1;
65 }
66 template<> int get(types::Bool SPARSE_CONST& d, int r, int c)
67 {
68     return d.get(r, c);
69 }
70 template<> bool get(types::SparseBool SPARSE_CONST& d, int r, int c)
71 {
72     return d.get(r, c);
73 }
74 template<> int get(types::SparseBool SPARSE_CONST& d, int r, int c)
75 {
76     return d.get(r, c);
77 }
78
79 template<> double get(types::Sparse SPARSE_CONST& s, int r, int c)
80 {
81     return s.getReal(r, c);
82 }
83 template<> std::complex<double> get(types::Sparse SPARSE_CONST& s, int r, int c)
84 {
85     return s.get(r, c);
86 }
87
88 template<> double get(types::Sparse::RealSparse_t SPARSE_CONST&s, int r, int c)
89 {
90     return s.coeff(r, c);
91 }
92 template<> std::complex<double> get(types::Sparse::RealSparse_t SPARSE_CONST&s, int r, int c)
93 {
94     return std::complex<double>(s.coeff(r, c), 0.);
95 }
96
97 template<> bool get(types::SparseBool::BoolSparse_t SPARSE_CONST& d, int r, int c)
98 {
99     return d.coeff(r, c);
100 }
101
102 template<> double get(types::Sparse::CplxSparse_t SPARSE_CONST&s, int r, int c)
103 {
104     return s.coeff(r, c).real();
105 }
106 template<> std::complex<double> get(types::Sparse::CplxSparse_t SPARSE_CONST&s, int r, int c)
107 {
108     return s.coeff(r, c);
109 }
110
111
112 /**
113    This free function overloads perform write access into a 2D container, using 0-based indices.
114    @param s the 2D structure used to fetch a value of type V.
115    @param r : the row (0 based)
116    @param c : the column (0 based)
117    @param v : the value of type V to set at row r and column c of structure s
118    @return : true iff everything went ok (should throw otherwise anyway).
119 */
120
121
122 template<typename S, typename V> bool set(S &, int, int, V)
123 {
124     return UndefinedAccessorForType<S>();
125 }
126
127 template<> bool set(types::Double & d, int r, int c, double v)
128 {
129     return d.set(r, c, v);
130 }
131 template<> bool set(types::Double & d, int r, int c, std::complex<double> v)
132 {
133     return d.set(r, c, v.real()) && d.setImg(r, c, v.imag());
134 }
135
136 template<> bool set(types::Sparse & s, int r, int c, double v)
137 {
138     return s.set(r, c, v);
139 }
140 template<> bool set(types::Sparse & s, int r, int c, std::complex<double> v)
141 {
142     return s.set(r, c, v);
143 }
144 template<> bool set(types::Bool & d, int r, int c, bool v)
145 {
146     return d.set(r, c, v);
147 }
148 template<> bool set(types::SparseBool & d, int r, int c, bool v)
149 {
150     return d.set(r, c, v);
151 }
152 template<> bool set(types::Bool & d, int r, int c, int v)
153 {
154     return d.set(r, c, v);
155 }
156 template<> bool set(types::SparseBool & d, int r, int c, int v)
157 {
158     return d.set(r, c, v != 0);
159 }
160
161 template<> bool set(types::Sparse::RealSparse_t& s, int r, int c, double v)
162 {
163     if (v != 0.)
164     {
165         s.insert(r, c) = v;
166     }
167     return true;
168 }
169
170 template<> bool set(types::Sparse::RealSparse_t& s, int r, int c, std::complex<double> v)
171 {
172     if ( v.real() != 0.)
173     {
174         s.insert(r, c) = v.real();
175     }
176     return  true;
177 }
178 // should we make this a compile error ?
179 template<> bool set(types::Sparse::CplxSparse_t& s, int r, int c, double v)
180 {
181     if (v != 0.)
182     {
183         s.insert(r, c) = std::complex<double>(v);
184     }
185     return true;
186 }
187
188 namespace
189 {
190 std::complex<double> const cplxZero(0., 0.);
191 }
192 template<> bool set(types::Sparse::CplxSparse_t& s, int r, int c, std::complex<double> v)
193 {
194     if (v != cplxZero)
195     {
196         s.insert(r, c) = v;
197     }
198     return true;
199 }
200
201 template<> bool set(types::SparseBool::BoolSparse_t& s, int r, int c, bool v)
202 {
203     if (v)
204     {
205         s.insert(r, c) = v;
206     }
207     return true;
208 }
209
210
211
212
213 template<typename S> inline int rows(S SPARSE_CONST&s)
214 {
215     return s.rows();
216 }
217 template<typename S> inline int cols(S SPARSE_CONST&s)
218 {
219     return s.cols();
220 }
221
222 template<> inline int rows(types::Double SPARSE_CONST&d)
223 {
224     return d.getRows();
225 }
226 template<> inline int cols(types::Double SPARSE_CONST&d)
227 {
228     return d.getCols();
229 }
230 template<> inline int rows(types::Sparse SPARSE_CONST&s)
231 {
232     return s.getRows();
233 }
234 template<> inline int cols(types::Sparse SPARSE_CONST&s)
235 {
236     return s.getCols();
237 }
238 template<> inline int rows(types::Bool SPARSE_CONST&s)
239 {
240     return s.getRows();
241 }
242 template<> inline int cols(types::Bool SPARSE_CONST&s)
243 {
244     return s.getCols();
245 }
246 template<> inline int rows(types::SparseBool SPARSE_CONST&s)
247 {
248     return s.getRows();
249 }
250 template<> inline int cols(types::SparseBool SPARSE_CONST&s)
251 {
252     return s.getCols();
253 }
254
255
256
257 /**
258   These free function overloads handle nb of rows size queries for 2D containers
259    wrapping the corresponding member function.
260    @param s : 2D structure to query
261    @return : nb of rows
262 */
263 template<typename S> inline int rows(S SPARSE_CONST&s);
264 template<> inline int rows(types::Double SPARSE_CONST&d);
265
266 /**
267   These free function overloads handle nb of cols size queries for 2D containers
268    wrapping the corresponding member function.
269    @param s : 2D structure to query
270    @return : nb of cols
271 */
272 template<typename S> inline int cols(S SPARSE_CONST&s);
273 template<> inline int cols(types::Double SPARSE_CONST&d);
274
275 /* this proxy struct provides read and write access (using set and get)
276    with the usual operators (operator*() and operator=() )*/
277 template<typename S, typename V> struct Accessor
278 {
279     /**
280        @param s_ : 2D structure to access
281        @param r_ : row to access
282        @param c_ ; column to access
283     */
284     Accessor(S& s_, int r_, int c_): s(s_), r(r_), c(c_) {}
285     /**
286        read accessor as a casting operator
287        @return : value of s at (r,c)
288      */
289     operator V() SPARSE_CONST
290     {
291         //        std::cerr<<"reading "<<get<S,V>(s, r, c)<<" @("<<r<<","<<c<<")\n";
292         return ::get<V>(s, r, c);
293     }
294     /**
295        write accessor as an assignment operator
296        @param v : value to set at (r,c) in s.
297     */
298     template<typename Sa, typename Va>
299     Accessor& operator=(Accessor<Sa, Va> const& a)
300     {
301         //        std::cerr<<"writing "<<( Va(const_cast<Accessor<Sa, Va>&>(a)))<<" @("<<r<<","<<c<<")\n";
302         //        Va tmp=const_cast<Accessor<Sa, Va>&>(a);
303         //        ::set<S,V>(s, r, c, tmp);
304         ::set<S, V>(s, r, c, Va(const_cast<Accessor<Sa, Va>&>(a)));
305         return *this;
306     }
307
308     Accessor& operator=(Accessor const& a)
309     {
310         //        std::cerr<<"writing "<<( V(const_cast<Accessor&>(a)))<<" @("<<r<<","<<c<<")\n";
311         ::set<S, V>(s, r, c, V(const_cast<Accessor&>(a)));
312         return *this;
313     }
314     Accessor& operator=(V const& v)
315     {
316         //        std::cerr<<"writing "<<v<<" @("<<r<<","<<c<<")\n";
317         ::set<S, V>(s, r, c, v);
318         return *this;
319     }
320 private:
321     S& s;
322     int r, c;
323 };
324
325 /* convenient typedef for pairs of (row, column) int values used as 2D coords */
326 typedef std::pair<int, int> Coords2D;
327 /* convenient typedef for iterator over pairs of (row, column) int values used as 2D coords */
328 typedef std::iterator<std::forward_iterator_tag, Coords2D > Coords2DIterator;
329 /**
330    Iterator over coords making a full row-wise traversal wrapping around when reaching
331    the end of the 2D container.
332  */
333 struct RowWiseFullIterator : Coords2DIterator
334 {
335     /**
336        @param cMax : size of the 2D structure
337      */
338     RowWiseFullIterator(Coords2D cMax): c(0, 0), cMax(cMax)
339     {
340     }
341     /**
342        @param cMax : size of the 2D structure
343        @param cInit : starting coords of the traversal.
344      */
345     RowWiseFullIterator(Coords2D cMax, Coords2D cInit): c(cInit), cMax(cMax)
346     {
347     }
348     /**
349        @param rm : nb of rows of the 2D structure
350        @param cm : nb of column of the 2D structure
351      */
352     RowWiseFullIterator(int rm, int cm): c(0, 0), cMax(rm, cm)
353     {
354     }
355     /**
356        @param rm : nb of rows of the 2D structure
357        @param cm : nb of column of the 2D structure
358        @param rInit : starting row of the traversal
359        @param cInit : starting column of the traversal
360      */
361     RowWiseFullIterator(int rm, int cm, int rInit, int cInit): c(rInit, cInit), cMax(rm, cm)
362     {
363     }
364     RowWiseFullIterator& operator++()
365     {
366         if (++c.first == cMax.first)
367         {
368             c.first = 0;
369             if (++c.second == cMax.second)
370             {
371                 /* wrap around */
372                 c.first = c.second = 0;
373             }
374         }
375         return *this;
376     }
377     RowWiseFullIterator operator++(int)
378     {
379         RowWiseFullIterator tmp(*this);
380         ++(*this);
381         return tmp;
382     }
383     std::pair<int, int> operator*() const
384     {
385         return c;
386     }
387 private:
388     Coords2D c;
389     Coords2D const cMax;
390 };
391
392 /**
393    Iterator over coords making a row-wise traversal of non zero elements of an Eigen Sparse Matrix
394  */
395 template<typename Sp>
396 struct RowWiseSparseIterator : Coords2DIterator
397 {
398     /**
399        @param sp: sparse matrix for non zero elements traversal
400      */
401     RowWiseSparseIterator(Sp const& sp): sp(sp), outerIdx(0), innerIt(sp, 0)
402     {
403     }
404     RowWiseSparseIterator& operator++()
405     {
406         ++innerIt;
407         if (!innerIt)
408         {
409             if (++outerIdx >= sp.outerSize())
410             {
411                 outerIdx = 0;
412             }
413             new (&innerIt) typename Sp::InnerIterator(sp, outerIdx);// innerIt= typename Sp::InnerIterator(sp, outerIdx) when Eigen will be fixed
414         }
415         return *this;
416     }
417     RowWiseSparseIterator operator++(int)
418     {
419         RowWiseFullIterator tmp(*this);
420         ++(*this);
421         return tmp;
422     }
423     std::pair<int, int> operator*() const
424     {
425         //        std::cerr<<"sparse it r="<<innerIt.row()<<" c="<<innerIt.col()<<std::endl;
426         return std::pair<int, int>(innerIt.row(), innerIt.col());
427     }
428 private:
429     Sp const& sp;
430     typename Eigen::internal::traits<Sp>::Index outerIdx;
431     typename Sp::InnerIterator innerIt;
432 };
433
434 /**
435    translate an iterator
436  */
437 template<typename C2DIter>
438 struct TranslatedIterator : Coords2DIterator
439 {
440     /**
441        @param C2DIter: translation as a vector of (rows, cols)
442        @param tr: translation as a vector of (rows, cols)
443      */
444     TranslatedIterator(C2DIter const& c2dIter, Coords2D tr): it(c2dIter), tr(tr)
445     {
446     }
447     TranslatedIterator& operator++()
448     {
449         ++it;
450         return *this;
451     }
452     TranslatedIterator operator++(int)
453     {
454         TranslatedIterator tmp(*this);
455         ++(*this);
456         return tmp;
457     }
458     std::pair<int, int> operator*() const
459     {
460         std::pair<int, int>res(*it);
461         res.first += tr.first;
462         res.second += tr.second;
463         //        std::cerr<<"translated it r="<< res.first<<" c="<<res.second<<std::endl;
464
465         return res;
466     }
467 private:
468     C2DIter it;
469     Coords2D const tr;
470 };
471
472 /**
473  * Template for iterator over 2D coords from an int*.
474  * Could handle wrap around with a length arg (i.e. to recycle values instead of raising
475  * "error 15 Submatrix incorrectly defined."
476  */
477 template<bool AsVector = false> struct Coords : Coords2DIterator
478 {
479     Coords(int SPARSE_CONST* coords, int unused = 0): coords(coords), unused(unused)
480     {
481     }
482
483     Coords& operator++()
484     {
485         coords += 2;
486         return *this;
487     }
488
489     Coords& operator++(int)
490     {
491         Coords tmp(*this);
492         ++(*this);
493         return tmp;
494     }
495
496     Coords2D operator*()const
497     {
498         return Coords2D(coords[0] - 1, coords[1] - 1);
499     }
500
501 private:
502     int const* coords;
503     int unused;
504 };
505 /**
506    explicit specialization for 2D from 1D int* sequences
507    (The 2D strcture is considered as a vector)
508  */
509 template<> struct Coords<true> : Coords2DIterator
510 {
511     Coords(int SPARSE_CONST* coords, int rMax): coords(coords), rMax(rMax)
512     {
513     }
514
515     Coords& operator++()
516     {
517         ++coords;
518         return *this;
519     }
520
521     Coords operator++(int)
522     {
523         Coords tmp(*this);
524         ++(*this);
525         return tmp;
526     }
527
528     Coords2D operator*()const
529     {
530         return Coords2D((coords[0] - 1) % rMax, (coords[0] - 1) / rMax);
531     }
532
533 private:
534     int const* coords;
535     int const rMax;
536 };
537 /* This 'iterator' class allows traverses the 2D containers, either
538 Rowwisefull traversal
539 or with 2D coords from another matrix
540 or with 1D coords from another vector (1x) matrix
541 to respect Double insert() API, we take int* and a bool
542 */
543 template<typename S, typename V, typename Iter>
544 struct MatrixIterator : std::iterator<std::forward_iterator_tag, V>
545 {
546     MatrixIterator(S& s_, Iter i_): s(s_), i(i_)
547     {
548     }
549     MatrixIterator& operator++()
550     {
551         ++i;
552         return *this;
553     }
554     MatrixIterator operator++(int)
555     {
556         MatrixIterator tmp(*this);
557         ++i;
558         return tmp;
559     }
560     Accessor<S, V> operator*()
561     {
562         return Accessor<S, V>(s, (*i).first, (*i).second);
563     }
564 private:
565     S& s;
566     Iter i;
567 };
568 template<typename V, typename S, typename Iter>
569 MatrixIterator<S, V, Iter> makeMatrixIterator(S& s, Iter i)
570 {
571     return MatrixIterator<S, V, Iter>(s, i);
572 }
573
574 template<typename S> struct IteratorFromVar;
575
576 template<typename S> IteratorFromVar<S> makeIteratorFromVar(S& s);
577
578 struct Adjacency
579 {
580     Adjacency(double const* x, double const*a): xadj(x), adjncy(a) {}
581     double const* xadj;
582     double const* adjncy;
583 };
584
585 template<typename In, typename Sz, typename Out>
586 Out mycopy_n(In i, Sz n, Out o)
587 {
588     for (; n; --n, ++i, ++o)
589     {
590         *o = *i;
591     }
592     return o;
593 }
594
595 template<typename T> std::size_t nonZeros(T SPARSE_CONST& t)
596 {
597     return t.getSize();
598 }
599 template<> std::size_t nonZeros(types::Sparse SPARSE_CONST& sp)
600 {
601     return sp.nonZeros();
602 }
603 template<typename Scalar, int Options, typename Index> std::size_t nonZeros(Eigen::SparseMatrix<Scalar, Options, Index> SPARSE_CONST& sp)
604 {
605     return sp.nonZeros();
606 }
607
608
609 /* Default for dense matrix Scilab matrix types
610  */
611 template<typename D> RowWiseFullIterator makeNonZerosIterator(D SPARSE_CONST& d)
612 {
613     return RowWiseFullIterator(d.getRows(), d.getCols());
614 }
615 template<typename Scalar, int Options, typename Index> RowWiseSparseIterator<Eigen::SparseMatrix<Scalar, Options, Index> > makeNonZerosIterator(Eigen::SparseMatrix<Scalar, Options, Index> SPARSE_CONST& sp)
616 {
617     return RowWiseSparseIterator<Eigen::SparseMatrix<Scalar, Options, Index> >(sp);
618 }
619 template<typename Iter> TranslatedIterator<Iter> makeTranslatedIterator(Iter const& it, Coords2D const& tr)
620 {
621     return TranslatedIterator<Iter>(it, tr);
622 }
623
624
625
626 template<typename S> struct IteratorFromVar { };
627
628 template<> struct IteratorFromVar<types::Double> : Coords2DIterator
629 {
630     IteratorFromVar(types::Double& d_): d(d_), r(0)
631     {
632         // check dimension ?
633     }
634
635     IteratorFromVar& operator++()
636     {
637         ++r;
638         return *this;
639     }
640     IteratorFromVar operator++(int)
641     {
642         IteratorFromVar tmp(*this);
643         ++r;
644         return tmp;
645     }
646     Coords2D operator*()
647     {
648         return std::pair<int, int>(static_cast<int>(d.getReal(r, 0) - 1), static_cast<int>(d.getReal(r, 1) - 1));
649     }
650 private:
651     types::Double& d;
652     int r;
653 };
654
655 /*
656   iterator from adjacency matrices :
657  */
658 template<> struct IteratorFromVar<Adjacency> : Coords2DIterator
659 {
660     IteratorFromVar(Adjacency& a): xadj(a.xadj), adjncy(a.adjncy), c(1), nb(1)
661     {
662         update();
663     }
664
665     IteratorFromVar& operator++()
666     {
667         ++nb;
668         update();
669         ++adjncy;
670         return *this;
671     }
672     IteratorFromVar operator++(int)
673     {
674         IteratorFromVar tmp(*this);
675         ++nb;
676         update();
677         ++adjncy;
678         return tmp;
679     }
680     std::pair<int, int> operator*()
681     {
682         return std::pair<int, int>(static_cast<int>(*adjncy) - 1, c - 1);
683     }
684 private:
685     void update()
686     {
687         for (; xadj[1] <= nb; ++c, ++xadj)
688         {
689         }
690     }
691     double const* xadj;
692     double const* adjncy;
693     int c;
694     std::size_t nb;
695 };
696
697 template<typename S> IteratorFromVar<S> makeIteratorFromVar(S& s)
698 {
699     return IteratorFromVar<S>(s);
700 }
701
702 #endif