1 /** @file gsGridIterator.h 2 3 @brief Provides iteration over integer or numeric points in a (hyper-)cube 4 5 This file is part of the G+Smo library. 6 7 This Source Code Form is subject to the terms of the Mozilla Public 8 License, v. 2.0. If a copy of the MPL was not distributed with this 9 file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 Author(s): A. Mantzaflaris 12 */ 13 14 #pragma once 15 16 namespace gismo 17 { 18 19 /** 20 \brief Specifies aliases describing the modes for gismo::gsGridIterator 21 22 \ingroup enums 23 */ 24 enum gsGridIteratorMode 25 { 26 CUBE = 0, ///< Cube mode iterates over all lattice points inside a cube 27 BDR = 1, ///< Boundary mode iterates over boundary lattice points only 28 VERTEX = 2, ///< Vertex mode iterates over cube vertices only 29 CWISE = 3 ///< Coordinate-wise mode iterates over a grid given by coordinate vectors 30 }; 31 32 // note: default arguments are found in gsForwardDeclarations.h 33 template<typename Z, int mode, short_t d, bool> class gsGridIterator 34 {using Z::GISMO_ERROR_gsGridIterator_has_invalid_template_arguments;}; 35 36 /** 37 \brief Iterator over the Cartesian product of integer points in a 38 tensor-product grid. 39 40 The iteration is done in lexicographic order. 41 42 - mode = 0 : iteration over [a, b) or [a, b] 43 gsGridIterator<T,CUBE> grid(a,b); 44 45 - mode = 1 : iteration over the boundry points of [a, b) or [a, b] 46 gsGridIterator<T,BDR> grid(a,b); 47 48 - mode = 2 : iteration over the vertices of [a, b) or [a, b] 49 gsGridIterator<T,VERTEX> grid(a,b); 50 51 52 The open or closed case is determined by a constructor flag. 53 A third argument, eg. gsGridIterator<T,CUBE,d> grid(a,b); 54 specifies fixed size dimension \a d 55 56 \note Iteration over the boundary including offsets is possible 57 using the free functions in gsUtils/gsCombinatorics.h 58 59 \tparam Z type of the integer coordinates of the index vector 60 61 \tparam d statically known dimension, or dynamic dimension if d = 62 -1 (default value) 63 64 \tparam mode 0: all integer points in [a,b], 1: all points in 65 [a,b), 2: vertices of cube [a,b] 66 67 \ingroup Tensor 68 */ 69 template<class Z, int mode, short_t d> 70 class gsGridIterator<Z,mode,d,true> 71 { 72 public: 73 typedef gsVector<Z,d> point; 74 75 public: 76 77 /** 78 \brief Empty constructor 79 */ gsGridIterator()80 gsGridIterator() 81 { 82 GISMO_STATIC_ASSERT(std::numeric_limits<Z>::is_integer,"The template parameter needs to be an integer type."); 83 GISMO_STATIC_ASSERT(mode > -1 && mode < 3, "The mode of gsGridIterator needs to be 0, 1 or 2."); 84 } 85 86 /** 87 \brief Constructor using lower and upper limits 88 89 \param a lower limit (vertex of a cube) 90 \param b upper limit (vertex of a cube) 91 \param open if true, the iteration is performed for the points in $[a_i,b_i)$ 92 */ 93 gsGridIterator(point const & a, point const & b, bool open = true) 94 { reset(a, b, open); } 95 96 /** 97 \brief Constructor using upper limit. The iteration starts from zero. 98 99 \param b upper limit (vertex of a cube) 100 \param open if true, the iteration is performed for the points in $[0,b_i)$ 101 */ 102 explicit gsGridIterator(point const & b, bool open = true) 103 { reset(point::Zero(b.size()), b, open); } 104 105 /** 106 \brief Constructor using lower and upper limits 107 108 \param ab Matrix with two columns, corresponding to lower and upper limits 109 \param open if true, the iteration is performed for the points in $[a_i,b_i)$ 110 */ 111 explicit gsGridIterator(gsMatrix<Z,d,2> const & ab, bool open = true) 112 { reset(ab.col(0), ab.col(1), open); } 113 114 /** 115 \brief Resets the iterator using new lower and upper limits 116 117 \param a lower limit (vertex of a cube) 118 \param b upper limit (vertex of a cube) 119 \param open if true, the iteration is performed for the points in $[a_i,b_i)$ 120 */ 121 inline void reset(point const & a, point const & b, bool open = true) 122 { 123 GISMO_ASSERT(a.rows() == b.rows(), "Invalid endpoint dimensions"); 124 m_low = m_cur = a; 125 if (open) m_upp = b.array() - 1; else m_upp = b; 126 m_valid = ( (m_low.array() <= m_upp.array()).all() ? true : false ); 127 } 128 129 /** 130 \brief Resets the iterator, so that a new iteration over the 131 points may start 132 */ reset()133 void reset() { reset(m_low,m_upp, false); } 134 135 // See http://eigen.tuxfamily.org/dox-devel/group__TopicStructHavingEigenMembers.html 136 //EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF( (sizeof(point)%16)==0 ); 137 138 public: 139 140 operator bool() const {return m_valid;} 141 142 const point & operator*() const {return m_cur;} 143 144 const point * operator->() const {return &m_cur;} 145 146 inline gsGridIterator & operator++() 147 { 148 // Multiple implementations according to mode 149 switch (mode) 150 { 151 case 0: // ----------- Iteration over [m_low, m_upp] 152 case 2: // iteration over vertices of the cube [m_low, m_upp] 153 for (index_t i = 0; i != m_cur.size(); ++i) 154 if ( m_cur[i] != m_upp[i] ) 155 { 156 if (0==mode) ++m_cur[i]; else m_cur[i] = m_upp[i]; 157 return *this; 158 } 159 else 160 m_cur[i] = m_low[i]; 161 m_valid = 0;//done 162 return *this; 163 164 case 1: // ----------- Iteration over boundary of [m_low, m_upp] 165 for (index_t i = 0; i != m_cur.size(); ++i) 166 { 167 if ( m_cur[i] != m_upp[i] ) 168 { 169 if ( m_cur[i] == m_low[i] && (i+1!=m_cur.size() || 1==m_cur.size()) ) 170 { 171 index_t c = i+1; 172 for (index_t j = c; j!=m_cur.size(); ++j) 173 if ( (m_cur[j] == m_low[j]) || 174 (m_cur[j] == m_upp[j]) ) 175 ++c; 176 177 if ( c==1 ) 178 m_cur[i] = m_upp[i]; 179 else 180 ++m_cur[i]; 181 } 182 else 183 m_cur[i]++; 184 185 m_cur.head(i) = m_low.head(i); 186 return *this; 187 } 188 } 189 /*fall through to default*/ 190 default: 191 m_valid = 0;//done 192 return *this; 193 } 194 } 195 196 /** 197 \brief Returns true if the \a i-th coordinate has minimal value 198 */ isFloor(int i)199 inline bool isFloor(int i) const { return m_cur[i] == m_low[i];} 200 201 /** 202 \brief Returns true if the \a i-th coordinate has maximal value 203 */ isCeil(int i)204 inline bool isCeil (int i) const 205 { return mode ? m_cur[i] == m_upp[i] : m_cur[i] + 1 == m_upp[i];} 206 207 /** 208 \brief Returns true if the current point lies on a boundary 209 */ isBoundary()210 bool isBoundary() const 211 { 212 return (m_cur.array() == m_low.array()).any() || 213 (m_cur.array() == m_upp.array()).any() ; 214 } 215 216 /** 217 \brief Returns the total number of points that are iterated 218 */ numPoints()219 index_t numPoints() const 220 { 221 switch (mode) 222 { 223 case 0: return ((m_upp-m_low).array() + 1).prod(); 224 case 1: return ((m_upp-m_low).array() + 1).prod() - ((m_upp-m_low).array()-1).prod(); 225 case 2: return ( 1 << m_low.rows() ); 226 default: GISMO_ERROR("gsGridIterator::numPoints"); 227 } 228 } 229 230 /** 231 \brief Returns the total number of points per coordinate which 232 are iterated 233 */ numPointsCwise()234 point numPointsCwise() const { return (m_upp-m_low).array() + 1;} 235 236 /** 237 \brief Returns the first point in the iteration 238 */ lower()239 const point & lower() const {return m_low;} 240 241 /** 242 \brief Returns the last point in the iteration 243 */ upper()244 const point & upper() const {return m_upp;} 245 246 /** 247 \brief Utility function which returns the vector of strides 248 249 \return An integer vector (stride vector). the entry \f$s_i\f$ 250 implies that the current \f$i\f$-th coordinate of the iterated 251 point will appear again after \f$s_i\f$ increments. Moreover, 252 the dot product of (*this - this->lower()) with the stride 253 vector results in the "flat index", i.e. the cardinal index of 254 the point in the iteration sequence. 255 256 If the iteration starts from the point zero 257 (this->lower()==zero), then the stride vector has the property 258 that that when we add it to *this we obtain the next point in 259 the iteration sequence. In this case the flat index is obtained 260 by taking the dot product with *this. 261 */ strides()262 point strides() const 263 { 264 point res(m_low.rows()); 265 res[0] = 1; 266 for (index_t i = 1; i != res.size(); ++i ) 267 res[i] = res[i-1] * ( m_upp[i-1] - m_low[i-1] + 1 ); 268 return res; 269 } 270 271 private: 272 273 point m_low, m_upp; ///< Iteration lower and upper limits 274 point m_cur; ///< Current point pointed at by the iterator 275 bool m_valid; ///< Indicates the state of the iterator 276 }; 277 278 279 /** 280 \brief Iterator over a Cartesian product of uniformly distributed 281 numeric points inside a (hyper-)cube. 282 283 The iteration is done in the natural lexicographic order. 284 285 - mode = 0 : iteration over uniform samples in [a, b] 286 gsGridIterator<T,CUBE> grid(a,b); 287 288 - mode = 1 : iteration over uniform samples on the boundary points [a, b] 289 gsGridIterator<T,BDR> grid(a,b); 290 291 - mode = 2 : iteration over the vertices [a, b] 292 gsGridIterator<T,VERTEX> grid(a,b); 293 294 A third argument, eg. gsGridIterator<T,CUBE,d> grid(a,b); 295 specifies fixed size dimension \a d. 296 297 \tparam T type of the numeric coordinates of the index vector 298 299 \tparam d statically known dimension, or dynamic dimension if d = 300 -1 (default value) 301 302 \tparam mode 0: all sample points in [a,b], 1: all points in 303 [a,b), 2: vertices of cube [a,b] 304 305 \ingroup Tensor 306 */ 307 template<class T, int mode, short_t d> 308 class gsGridIterator<T,mode,d,false> 309 { 310 public: 311 312 typedef gsVector<T,d> point; 313 314 typedef gsGridIterator<index_t, mode, d> integer_iterator; 315 316 typedef typename integer_iterator::point point_index; 317 public: 318 319 /** 320 \brief Constructor using lower and upper limits 321 322 Uniformly sampled points will be generated 323 324 \param a lower limit (vertex of a cube) 325 \param b upper limit (vertex of a cube) 326 \param np number of sample points per coordinate 327 */ gsGridIterator(point const & a,point const & b,point_index const & np)328 gsGridIterator(point const & a, 329 point const & b, 330 point_index const & np) 331 : m_iter(np, 1) 332 { 333 reset(a, b); 334 GISMO_STATIC_ASSERT(mode > -1 && mode < 3, "The mode of gsGridIterator needs to be 0, 1 or 2."); 335 } 336 337 /** 338 \brief Constructor using lower and upper limits 339 340 Uniformly sampled points will be generated 341 342 \param ab Matrix with two columns, corresponding to lower and upper limits 343 \param np number of sample points per coordinate 344 */ gsGridIterator(gsMatrix<T,d,2> const & ab,point_index const & np)345 gsGridIterator(gsMatrix<T,d,2> const & ab, 346 point_index const & np) 347 : m_iter(np, 1) 348 { 349 reset(ab.col(0), ab.col(1)); 350 } 351 352 /** 353 \brief Constructor using lower and upper limits 354 355 Uniformly sampled points will be generated. The number of the 356 points is approximately \a numPoints 357 358 \param ab Matrix with two columns, corresponding to lower and upper limits 359 \param numPoints the number sample points to be used (in total, approximately) 360 */ gsGridIterator(gsMatrix<T,d,2> const & ab,unsigned numPoints)361 gsGridIterator(gsMatrix<T,d,2> const & ab, unsigned numPoints) 362 : m_low(ab.col(0)), m_upp(ab.col(1)) 363 { 364 // deduce the number of points per direction for an approximately uniform grid 365 const gsVector<double,d> L = (ab.col(1) - ab.col(0)).template cast<double>(); 366 const double h = std::pow(L.prod()/numPoints, 1.0 / m_low.rows()); 367 const point_index npts = (L/h).unaryExpr((double(*)(double))std::ceil).template cast<index_t>(); 368 m_iter = integer_iterator(npts, 1); 369 reset(ab.col(0), ab.col(1)); 370 } 371 372 /** 373 \brief Resets the iterator, so that a new iteration over the 374 points may start 375 */ reset()376 void reset() { m_cur = m_low; m_iter.reset(); } 377 378 /** 379 \brief Resets the iterator using new lower and upper limits 380 381 \param a lower limit (vertex of a cube) 382 \param b upper limit (vertex of a cube) 383 */ reset(point const & a,point const & b)384 void reset(point const & a, 385 point const & b) 386 { 387 m_cur = m_low = a; 388 m_upp = b; 389 m_step = (b-a).array() / (m_iter.numPointsCwise().array() - 1) 390 .matrix().cwiseMax(1).template cast<T>().array() ; 391 m_iter.reset(); 392 } 393 toMatrix()394 gsMatrix<T> toMatrix() const 395 { 396 gsMatrix<T> res(m_cur.rows(), numPoints()); 397 integer_iterator iter = m_iter; 398 iter.reset(); 399 400 for(index_t c = 0; iter; ++iter, ++c) 401 update(*iter,res.col(c).data()); 402 403 return res; 404 } 405 406 // See http://eigen.tuxfamily.org/dox-devel/group__TopicStructHavingEigenMembers.html 407 //EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF( (sizeof(point)%16)==0 ); 408 409 public: 410 411 operator bool() const {return m_iter;} 412 413 const gsMatrix<T> & operator*() const {return m_cur;} 414 415 const gsMatrix<T> * operator->() const {return &m_cur;} 416 417 inline gsGridIterator & operator++() 418 { 419 if ( ++m_iter ) update(*m_iter, m_cur.data()); 420 return *this; 421 } 422 423 /** 424 \brief Returns true if the \a i-th coordinate has minimal value 425 */ isFloor(int i)426 inline bool isFloor(int i) const { return m_iter.isFloor(i);} 427 428 /** 429 \brief Returns true if the \a i-th coordinate has maximal value 430 */ isCeil(int i)431 inline bool isCeil (int i) const { return m_iter.isCeil(i);} 432 433 /** 434 \brief Returns true if the current point lies on a boundary 435 */ isBoundary()436 inline bool isBoundary() const { return m_iter.isBoundary();} 437 438 /** 439 \brief Returns the total number of points that are iterated 440 */ numPoints()441 index_t numPoints() const {return m_iter.numPoints();} 442 443 /** 444 \brief Returns the total number of points per coordinate which 445 are iterated 446 */ numPointsCwise()447 point_index numPointsCwise() const {return m_iter.numPointsCwise();} 448 449 /** 450 \brief Returns the first point in the iteration 451 */ lower()452 const point & lower() const {return m_low;} 453 454 /** 455 \brief Returns the last point in the iteration 456 */ upper()457 const point & upper() const {return m_upp;} 458 459 /** 460 \brief Returns the tensor index of the current point 461 */ tensorIndex()462 const point_index & tensorIndex() const { return *m_iter;} 463 464 /** 465 \brief Returns the \a i-th index of the current point 466 */ index(const index_t i)467 index_t index(const index_t i) const { return m_iter->at(i);} 468 469 /** 470 \brief Returns the coordinate-wise stepping between the samples 471 */ step()472 const point & step() const {return m_step;} 473 474 /** 475 \brief Returns a reference to the underlying integer lattice 476 iterator 477 */ index_iterator()478 const integer_iterator & index_iterator() const { return m_iter;} 479 480 /** 481 \brief Returns the point corresponding to tensor index \a ti 482 483 \note Unlikely to std iterators, the position is not 484 counted from the current position of \a this, but from the 485 starting point of the iteration 486 487 \param ti a valid tensor index of a point in the iteration sequence 488 */ 489 inline const gsMatrix<T> operator [] (const point_index & ti) const 490 { 491 gsMatrix<T> res(m_low.rows(),1); 492 update(ti, res.data()); 493 return res; 494 } 495 496 private: 497 498 // Update the point \a pt to the position \a ti update(const point_index & ti,T * pt)499 inline void update(const point_index & ti, T * pt) const 500 { 501 for (index_t i = 0; i != m_low.size(); ++i) 502 if ( ti[i] == m_iter.lower()[i] ) 503 *(pt++) = m_low[i]; // avoid numerical error at first val 504 else if ( ti[i] == m_iter.upper()[i] ) 505 *(pt++) = m_upp[i]; // avoid numerical error at last val 506 else 507 *(pt++) = m_low[i] + ti[i] * m_step[i]; 508 } 509 510 private: 511 512 point m_low, m_upp, m_step; ///< Iteration lower and upper limits and stepsize 513 514 integer_iterator m_iter; ///< Underlying integer lattice iterator 515 gsMatrix<T> m_cur; ///< Current point pointed at by the iterator 516 }; 517 518 519 /** 520 \brief Iterator over a Cartesian product of points, which is 521 given by coordinate-wise point sets 522 523 The iteration is done in lexicographic order. Usage: 524 525 gsGridIterator<T,CWISE> grid(a,b); 526 527 or 528 529 gsGridIterator<T,CWISE,d> grid(a,b); 530 531 \tparam T type of the numeric coordinates of the index vector 532 533 \tparam d statically known dimension, or dynamic dimension if d = 534 -1 (default value) 535 536 \ingroup Tensor 537 */ 538 template<class T, short_t d> // mode == 3 == CWISE 539 class gsGridIterator<T,CWISE,d,false> 540 { 541 public: 542 543 typedef gsGridIterator<index_t, 0, d> integer_iterator; 544 545 typedef typename integer_iterator::point point_index; 546 547 typedef gsVector<const T *,d, 2> CwiseData; //non-aligned array 548 549 public: 550 551 /** 552 \brief Constructor using references to the coordinate vectors 553 554 \param cwise A container of matrices or vectors, each 555 containing the sample points in the respective coordinate 556 */ 557 template<class CwiseContainer> gsGridIterator(const CwiseContainer & cwise)558 explicit gsGridIterator(const CwiseContainer & cwise) 559 : m_cwise(cwise.size()), m_cur(m_cwise.size(),1) 560 { 561 point_index npts(m_cwise.size()); 562 for (index_t i = 0; i != npts.size(); ++i) 563 { 564 m_cwise[i] = cwise[i].data(); 565 npts[i] = cwise[i].size() - 1; 566 GISMO_ASSERT(cwise[i].cols()==1 || cwise[i].rows()==1, "Invalid input"); 567 } 568 m_iter = integer_iterator(npts, 0); 569 //if (1==cwise.front().cols()) 570 // m_cur.derived().resize(1, npts.size()); 571 //else 572 //m_cur.derived().resize(npts.size(), 1); 573 update(*m_iter, m_cur.data()); 574 } 575 576 template<class Matrix_t> gsGridIterator(const std::vector<Matrix_t * > & cwise)577 explicit gsGridIterator(const std::vector<Matrix_t*> & cwise) 578 : m_cwise(cwise.size()), m_cur(m_cwise.size(),1) 579 { 580 point_index npts(m_cwise.size()); 581 for (index_t i = 0; i != npts.size(); ++i) 582 { 583 m_cwise[i] = cwise[i]->data(); 584 npts[i] = cwise[i]->size() - 1; 585 GISMO_ASSERT(cwise[i]->cols()==1 || cwise[i]->rows()==1, "Invalid input"); 586 } 587 m_iter = integer_iterator(npts, 0); 588 update(*m_iter, m_cur.data()); 589 } 590 591 /** 592 \brief Resets the iterator, so that a new iteration over the 593 points may start 594 */ reset()595 void reset() { m_iter.reset(); update(*m_iter, m_cur.data());} 596 597 //void restart() { m_iter.reset(); update(*m_iter, m_cur.data());} 598 599 public: 600 601 operator bool() const {return m_iter;} 602 603 const gsMatrix<T> & operator*() const {return m_cur;} 604 605 const gsMatrix<T> * operator->() const {return &m_cur;} 606 607 inline gsGridIterator & operator++() 608 { 609 if (++m_iter) update(*m_iter, m_cur.data()); 610 return *this; 611 } 612 613 /** 614 \brief Returns true if the \a i-th coordinate has minimal value 615 */ isFloor(int i)616 inline bool isFloor(int i) const { return m_iter.isFloor(i);} 617 618 /** 619 \brief Returns true if the \a i-th coordinate has maximal value 620 */ isCeil(int i)621 inline bool isCeil (int i) const { return m_iter.isCeil(i);} 622 623 /** 624 \brief Returns true if the current point lies on a boundary 625 */ isBoundary()626 inline bool isBoundary() const { return m_iter.isBoundary();} 627 628 /** 629 \brief Returns the total number of points that are iterated 630 */ numPoints()631 index_t numPoints() const {return m_iter.numPoints();} 632 633 /** 634 \brief Returns the total number of points per coordinate which 635 are iterated 636 */ numPointsCwise()637 point_index numPointsCwise() const {return m_iter.numPointsCwise();} 638 639 /** 640 \brief Returns the tensor index of the current point 641 */ tensorIndex()642 const point_index & tensorIndex() const { return *m_iter;} 643 644 /** 645 \brief Returns the \a i-th index of the current point 646 */ index(const index_t i)647 index_t index(const index_t i) const { return m_iter->at(i);} 648 649 /** 650 \brief Returns a reference to the underlying integer lattice 651 iterator 652 */ index_iterator()653 const integer_iterator & index_iterator() const { return m_iter;} 654 655 /** 656 \brief Returns the point corresponding to tensor index \a ti 657 658 \note Unlikely to std iterators, the position is not 659 counted from the current position of \a this, but from the 660 starting point of the iteration 661 662 \param ti a valid tensor index of a point in the iteration sequence 663 */ 664 inline const gsMatrix<T> operator [] (const point_index & ti) const 665 { 666 gsMatrix<T> res(m_cur.rows(),1); 667 update(ti, res.data()); 668 return res; 669 } 670 671 private: 672 673 // Update the point \a pt to the position \a ti update(const point_index & ti,T * pt)674 inline void update(const point_index & ti, T * pt) 675 { 676 for (index_t i = 0; i != m_cur.rows(); ++i) 677 *(pt++) = m_cwise[i][ti[i]]; 678 } 679 680 private: 681 682 CwiseData m_cwise; ///< List of coordinate-wise values 683 684 integer_iterator m_iter; ///< Underlying integer lattice iterator 685 gsMatrix<T> m_cur; ///< Current point pointed at by the iterator 686 }; 687 688 689 } // namespace gismo 690