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