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