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