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