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