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