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 /*
 161  * TODO report possible bug in Eigen when inserting 0. invalidates Eigen::InnerIterator
 162  */
 163 template<> bool set(types::Sparse::RealSparse_t& 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(types::Sparse::RealSparse_t& 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(types::Sparse::CplxSparse_t& 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(types::Sparse::CplxSparse_t& 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(types::SparseBool::BoolSparse_t& 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