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)
  48 {
  49     return UndefinedAccessorForType<S>();
  50 }
  51
  52 template<> double get(types::Double CONST& d, int r, int c)
  53 {
  54     return d.getReal(r, c);
  55 }
  56 template<> std::complex<double> get(types::Double 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 CONST& d, int r, int c)
  62 {
  63     return d.get(r, c) == 1;
  64 }
  65 template<> int get(types::Bool CONST& d, int r, int c)
  66 {
  67     return d.get(r, c);
  68 }
  69 template<> bool get(types::SparseBool CONST& d, int r, int c)
  70 {
  71     return d.get(r, c);
  72 }
  73 template<> int get(types::SparseBool CONST& d, int r, int c)
  74 {
  75     return d.get(r, c);
  76 }
  77
  78 template<> double get(types::Sparse CONST& s, int r, int c)
  79 {
  80     return s.getReal(r, c);
  81 }
  82 template<> std::complex<double> get(types::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> 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> 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> 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> 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> 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 CONST&s)
 216 {
 217     return s.rows();
 218 }
 219 template<typename S> int cols(S CONST&s)
 220 {
 221     return s.cols();
 222 }
 223
 224 template<> int rows(types::Double CONST&d)
 225 {
 226     return d.getRows();
 227 }
 228 template<> int cols(types::Double CONST&d)
 229 {
 230     return d.getCols();
 231 }
 232 template<> int rows(types::Sparse CONST&s)
 233 {
 234     return s.getRows();
 235 }
 236 template<> int cols(types::Sparse CONST&s)
 237 {
 238     return s.getCols();
 239 }
 240 template<> int rows(types::Bool CONST&s)
 241 {
 242     return s.getRows();
 243 }
 244 template<> int cols(types::Bool CONST&s)
 245 {
 246     return s.getCols();
 247 }
 248 template<> int rows(types::SparseBool CONST&s)
 249 {
 250     return s.getRows();
 251 }
 252 template<> int cols(types::SparseBool 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 CONST&s);
 266 template<> inline int rows(types::Double 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 CONST&s);
 275 template<> inline int cols(types::Double 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() 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 CONST* coords, int unused = 0): coords(coords)
 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 };
 506 /**
 507    explicit specialization for 2D from 1D int* sequences
 508    (The 2D strcture is considered as a vector)
 509  */
 510 template<> struct Coords<true> : Coords2DIterator
 511 {
 512     Coords(int CONST* coords, int rMax): coords(coords), rMax(rMax)
 513     {
 514     }
 515
 516     Coords& operator++()
 517     {
 518         ++coords;
 519         return *this;
 520     }
 521
 522     Coords operator++(int)
 523     {
 524         Coords tmp(*this);
 525         ++(*this);
 526         return tmp;
 527     }
 528
 529     Coords2D operator*()const
 530     {
 531         return Coords2D((coords[0] - 1) % rMax, (coords[0] - 1) / rMax);
 532     }
 533
 534 private:
 535     int const* coords;
 536     int const rMax;
 537 };
 538 /* This 'iterator' class allows traverses the 2D containers, either
 539 Rowwisefull traversal
 540 or with 2D coords from another matrix
 541 or with 1D coords from another vector (1x) matrix
 542 to respect Double insert() API, we take int* and a bool
 543 */
 544 template<typename S, typename V, typename Iter>
 545 struct MatrixIterator : std::iterator<std::forward_iterator_tag, V>
 546 {
 547     MatrixIterator(S& s_, Iter i_): s(s_), i(i_)
 548     {
 549     }
 550     MatrixIterator& operator++()
 551     {
 552         ++i;
 553         return *this;
 554     }
 555     MatrixIterator operator++(int)
 556     {
 557         MatrixIterator tmp(*this);
 558         ++i;
 559         return tmp;
 560     }
 561     Accessor<S, V> operator*()
 562     {
 563         return Accessor<S, V>(s, (*i).first, (*i).second);
 564     }
 565 private:
 566     S& s;
 567     Iter i;
 568 };
 569 template<typename V, typename S, typename Iter>
 570 MatrixIterator<S, V, Iter> makeMatrixIterator(S& s, Iter i)
 571 {
 572     return MatrixIterator<S, V, Iter>(s, i);
 573 }
 574
 575 template<typename S> struct IteratorFromVar;
 576
 577 template<typename S> IteratorFromVar<S> makeIteratorFromVar(S& s);
 578
 579 struct Adjacency
 580 {
 581     Adjacency(double const* x, double const*a): xadj(x), adjncy(a) {}
 582     double const* xadj;
 583     double const* adjncy;
 584 };
 585
 586 template<typename In, typename Sz, typename Out>
 587 Out mycopy_n(In i, Sz n, Out o)
 588 {
 589     for (; n; --n, ++i, ++o)
 590     {
 591         *o = *i;
 592     }
 593     return o;
 594 }
 595
 596 template<typename T> std::size_t nonZeros(T CONST& t)
 597 {
 598     return t.getSize();
 599 }
 600 template<> std::size_t nonZeros(types::Sparse CONST& sp)
 601 {
 602     return sp.nonZeros();
 603 }
 604 template<typename Scalar, int Options, typename Index> std::size_t nonZeros(Eigen::SparseMatrix<Scalar, Options, Index> CONST& sp)
 605 {
 606     return sp.nonZeros();
 607 }
 608
 609
 610 /* Default for dense matrix Scilab matrix types
 611  */
 612 template<typename D> RowWiseFullIterator makeNonZerosIterator(D CONST& d)
 613 {
 614     return RowWiseFullIterator(d.getRows(), d.getCols());
 615 }
 616 template<typename Scalar, int Options, typename Index> RowWiseSparseIterator<Eigen::SparseMatrix<Scalar, Options, Index> > makeNonZerosIterator(Eigen::SparseMatrix<Scalar, Options, Index> CONST& sp)
 617 {
 618     return RowWiseSparseIterator<Eigen::SparseMatrix<Scalar, Options, Index> >(sp);
 619 }
 620 template<typename Iter> TranslatedIterator<Iter> makeTranslatedIterator(Iter const& it, Coords2D const& tr)
 621 {
 622     return TranslatedIterator<Iter>(it, tr);
 623 }
 624
 625
 626
 627 template<typename S> struct IteratorFromVar { };
 628
 629 template<> struct IteratorFromVar<types::Double> : Coords2DIterator
 630 {
 631     IteratorFromVar(types::Double& d_): d(d_), r(0)
 632     {
 633         // check dimension ?
 634     }
 635
 636     IteratorFromVar& operator++()
 637     {
 638         ++r;
 639         return *this;
 640     }
 641     IteratorFromVar operator++(int)
 642     {
 643         IteratorFromVar tmp(*this);
 644         ++r;
 645         return tmp;
 646     }
 647     Coords2D operator*()
 648     {
 649         return std::pair<int, int>(static_cast<int>(d.getReal(r, 0) - 1), static_cast<int>(d.getReal(r, 1) - 1));
 650     }
 651 private:
 652     types::Double& d;
 653     int r;
 654 };
 655
 656 /*
 657   iterator from adjacency matrices :
 658  */
 659 template<> struct IteratorFromVar<Adjacency> : Coords2DIterator
 660 {
 661     IteratorFromVar(Adjacency& a): xadj(a.xadj), adjncy(a.adjncy), c(1), nb(1)
 662     {
 663         update();
 664     }
 665
 666     IteratorFromVar& operator++()
 667     {
 668         ++nb;
 669         update();
 670         ++adjncy;
 671         return *this;
 672     }
 673     IteratorFromVar operator++(int)
 674     {
 675         IteratorFromVar tmp(*this);
 676         ++nb;
 677         update();
 678         ++adjncy;
 679         return tmp;
 680     }
 681     std::pair<int, int> operator*()
 682     {
 683         return std::pair<int, int>(static_cast<int>(*adjncy) - 1, c - 1);
 684     }
 685 private:
 686     void update()
 687     {
 688         for (; xadj[1] <= nb; ++c, ++xadj)
 689         {
 690         }
 691     }
 692     double const* xadj;
 693     double const* adjncy;
 694     int c;
 695     std::size_t nb;
 696 };
 697
 698 template<typename S> IteratorFromVar<S> makeIteratorFromVar(S& s)
 699 {
 700     return IteratorFromVar<S>(s);
 701 }
 702
 703 #endif