MultiIndex.hpp

00001 #ifndef __Fiber_MultiIndex_HPP
00002 #define __Fiber_MultiIndex_HPP
00003 
00004 #include <eagle/vecmath.hpp>
00005 #include <eagle/Assignment.hpp>
00006 #include <eagle/Operator.hpp>
00007 #include <eagle/FixedArray.hpp>
00008 #include "Iterator.hpp"
00009 
00010 namespace Fiber
00011 {
00012 
00037 struct  CreateFromLinear {};
00038 
00039 struct  Power2Alignment {};
00040 
00052 template <int Dims>
00053         class MultiIndex : public MultiIndex<Dims-1>
00054 {
00055         typedef MultiIndex<Dims-1>  base_t;
00056 
00057         index_t idx;
00058 
00059 public:
00060         typedef index_t value_type;
00061 
00063         MultiIndex(const MultiIndex&M, const MultiIndex&D, const Add&) throw()
00064         : base_t(M,D, Add() )
00065         , idx(M.idx + D.maxidx())
00066         {}
00067 
00069         MultiIndex(const MultiIndex&M, const MultiIndex&D, const Sub&) throw()
00070         : base_t(M,D, Sub() )
00071         , idx(M.idx - D.maxidx())
00072         {}
00073 
00075         MultiIndex(const MultiIndex&M, const MultiIndex&D, const Mult&) throw()
00076         : base_t(M,D, Mult() )
00077         , idx(M.idx * D.maxidx())
00078         {}
00079 
00081         MultiIndex(const MultiIndex&M, const MultiIndex&D, const Div&)
00082         : base_t(M,D, Div() )
00083         , idx(M.idx / D.maxidx())
00084         {}
00085 
00086         MultiIndex(unsigned int bits, const ::Eagle::BinaryAnd&) throw()
00087         : base_t(bits, ::Eagle::BinaryAnd() )
00088         , idx( bits & (1<<(Dims-1)) ? 1 : 0 )
00089         {}
00090 
00091         MultiIndex(const MultiIndex&M, const MultiIndex&D, const ::Eagle::BinaryAnd&) throw()
00092         : base_t(M,D, ::Eagle::BinaryAnd() )
00093         , idx( M.idx > D.maxidx() ? M.idx - D.maxidx() : D.maxidx() - M.idx )
00094         {}
00095 
00096         MultiIndex(const MultiIndex&M, const Power2Alignment&) throw()
00097         : base_t(M, Power2Alignment() )
00098         , idx( RoundUpToNextHighestPowerOf2( M.idx) )
00099         {}
00100 
00101 public:
00106         friend MultiIndex distance(const MultiIndex&M0, const MultiIndex&M1)
00107         {
00108                 return MultiIndex(M0,M1, ::Eagle::BinaryAnd() );
00109         }
00110 
00111         friend MultiIndex bitalign(const MultiIndex&M) throw()
00112         {
00113                 return MultiIndex(M, Power2Alignment() );
00114         }
00115 
00135 static  MultiIndex BitIndex(unsigned int bits) throw()
00136         {
00137                 return MultiIndex(bits, ::Eagle::BinaryAnd() );
00138         }
00139 
00150 static  MultiIndex Axis(unsigned int orientation) throw()
00151         {
00152                 return BitIndex( 1 << orientation );
00153         }
00154 
00170         index_t BitIndex() const throw()
00171         {
00172                 return linear( MultiIndex(2) );
00173         }
00174 
00179 static  int     log2(index_t N) throw()
00180         {
00181                 if (N & (1<<(Dims-1) ) )
00182                     return Dims-1;
00183 
00184                 return base_t::log2(N);
00185         }
00186 
00206         index_t Orientation() const throw()
00207         {
00208                 return log2( BitIndex() );
00209         }       
00210 
00211         enum
00212         {
00214                 dims = Dims,
00215 
00217                 SIZE = Dims
00218         };
00219 
00229         int getSignedOrientation(const MultiIndex&B) const
00230         {
00231                 if (idx<B.idx) return +Dims;
00232                 if (idx>B.idx) return -Dims;
00233                 return subidx().getSignedOrientation(B.subidx() );
00234         }
00235 
00237         MultiIndex() throw()
00238         : idx(0)
00239         {}
00240 
00242         explicit MultiIndex(const index_t&I) throw()
00243         : base_t(I), idx(I)
00244         {}
00245 
00247         MultiIndex(const MultiIndex<Dims-1>&Midx, const index_t&I) throw()
00248         : base_t(Midx), idx(I)
00249         {}
00250 
00256         MultiIndex(const MultiIndex<Dims-1>&m, const index_t&Slice, int SliceDirection)
00257         {
00258                 for(int i=0, k=0; i<Dims; i++)
00259                 {
00260                         if (i==SliceDirection)
00261                                 (*this)[i] = Slice; 
00262                         else
00263                         {
00264                                 (*this)[i] = m[k]; 
00265                                 k++;
00266                         }
00267                 }
00268         }
00269 
00290         MultiIndex(const index_t&lin_idx, const MultiIndex&Dimens, const CreateFromLinear&) throw()
00291         : base_t( lin_idx, Dimens, CreateFromLinear() )
00292         , idx( (lin_idx / Dimens.subidx().size()) % Dimens.maxidx() )
00293         {}
00294 
00296         const MultiIndex<Dims-1>&subidx() const throw()
00297         {       return *this;   }
00298 
00300         MultiIndex<Dims-1>&subidx() throw()
00301         {       return *this;   }
00302 
00308         const index_t&operator[](int i) const
00309         {       return (i>=Dims-1) ? idx : subidx()[i];         }
00310 
00314         index_t&operator[](int i) throw()
00315         {       return (i>=Dims-1) ? idx : subidx()[i];         }
00316 
00317 
00319         Eagle::Assignment<MultiIndex,0> operator=(const index_t&i)
00320         {
00321                 return Eagle::Assignment<MultiIndex, 0>(*this, i);
00322         }
00323 
00324         const index_t&maxidx() const throw()
00325         {       return idx; }
00326 
00327         index_t&maxidx() throw()
00328         {       return idx; }
00329 
00331         index_t size() const throw()
00332         {
00333                 return idx*subidx().size();     
00334         }
00335 
00337         bool operator!=(const MultiIndex&D) const throw()
00338         {
00339                 return idx != D.maxidx() || subidx() != D.subidx() ;
00340         }
00341 
00343         bool operator==(const MultiIndex&M) const throw()
00344         {
00345                 return idx == M.maxidx() && subidx() == M.subidx() ;
00346         }
00347 
00349         MultiIndex&operator+=(const MultiIndex&D) throw()
00350         {
00351                 idx += D.maxidx(); 
00352                 subidx() += D.subidx(); 
00353                 return *this;
00354         }
00355 
00357         MultiIndex&operator-=(const MultiIndex&D) throw()
00358         {
00359                 idx -= D.maxidx(); 
00360                 subidx() -= D.subidx();
00361                 return *this;
00362         }
00363 
00364 /*
00365         friend MultiIndex clamp(const MultiIndex&M, const MultiIndex&D)
00366         {
00367         MultiIndex retval; 
00368                 retval.idx      = clamp(M.idx     , D.idx);
00369                 retval.subidx() = clamp(M.subidx(), D.subidx()); 
00370                 return retval;
00371         }
00372 
00373         friend MultiIndex clamp_m1(const MultiIndex&M, const MultiIndex&D)
00374         {
00375         MultiIndex retval;
00376                 retval.idx      = clamp_m1(M.idx     , D.idx-1);
00377                 retval.subidx() = clamp_m1(M.subidx(), D.subidx());
00378                 return retval;
00379         }
00380 */
00381 
00387         bool    inc(const MultiIndex&Dimens) throw()
00388         {
00389                 idx ++;
00390                 if ( idx >= Dimens.maxidx() )
00391                 {
00392                 bool    inc = subidx().inc(Dimens);
00393                         if (inc)
00394                                 idx = index_t(0); 
00395                         return inc;
00396                 } 
00397                 return true;
00398         }
00399 
00406         bool    inc(const MultiIndex&Dimens, const MultiIndex&Increment) throw()
00407         {
00408                 idx += Increment.maxidx();
00409                 if ( idx >= Dimens.maxidx() )
00410                 {
00411                 bool    inc = subidx().inc( Dimens, Increment );
00412                         if (inc)
00413                                 idx = index_t(0); 
00414                         return inc;
00415                 } 
00416                 return true;
00417         }
00418 
00422         index_t linear(const MultiIndex&Dimens) const throw()
00423         {       
00424 //              return idx + Dimens.maxidx() * subidx().linear( Dimens.subidx() ); 
00425 //              return idx * Dimens.maxidx() + subidx().linear( Dimens.subidx() ); 
00426 //              return idx*subidx().size(); 
00427 
00428                 assert(idx >= 0 && idx < Dimens.maxidx() );
00429 
00430                 /*
00431                  Pretty sure, this formula can be optimized.
00432                  The size() call on the subdimension recursively multiplies
00433                  its elements, but this could be moved into the summation/
00434                  multiplication recursive call.
00435                 */
00436         index_t L = Dimens.subidx().size();
00437                 return idx * L + subidx().linear( Dimens.subidx() ); 
00438         }
00439 
00440 
00441 //               MSVC                           GCC             G++                   SGI C++
00442 #if     defined(_IOSTREAM_) || defined(_GLIBCXX_OSTREAM) || _CPP_IOSTREAM || defined(__SGI_STL_IOSTREAM)
00443 
00444         std::ostream&print(std::ostream&os, char sep) const
00445         {       
00446                 base_t::print(os, sep); 
00447                 return os << sep << idx;
00448         }
00449 
00451         friend std::ostream&operator<<(std::ostream&os, const MultiIndex&MI)
00452         {
00453                 return MI.print(os, 'x');
00454         }
00455 #endif
00456 
00458         bool    isWithin(const MultiIndex&Range) const throw()
00459         {
00460         // return maxidx()>=0 && maxidx() < Range.maxidx() &&
00461                 return  maxidx() < Range.maxidx() &&
00462                         subidx().isWithin( Range );
00463         }
00464 
00466         bool    operator<(const MultiIndex&Range) const throw()
00467         {
00468                 return isWithin( Range );
00469         }
00470 
00472         bool    operator>(const MultiIndex&Range) const throw()
00473         {
00474                 return !isWithin( Range );
00475         }
00476 
00478         MultiIndex operator+(const MultiIndex&D) const throw()
00479         {       return MultiIndex(*this, D, Add() );    }
00480 
00482         MultiIndex operator-(const MultiIndex&D) const throw()
00483         {       return MultiIndex(*this, D, Sub() );    }
00484 
00486         MultiIndex operator*(const MultiIndex&D) const throw()
00487         {       return MultiIndex(*this, D, Mult() );   }
00488 
00490         MultiIndex operator/(const MultiIndex&D) const
00491         {       return MultiIndex(*this, D, Div() );    }
00492 
00493 
00503         MultiIndex operator+(int i) const
00504         {
00505         const MultiIndex&M = *this;
00506                 if (i==0) return M;
00507 
00508                 if (i<0)
00509                         return M - MultiIndex::Axis(-i+1); 
00510                 else
00511                         return M + MultiIndex::Axis( i-1);
00512         }
00513 
00519         friend MultiIndex operator-(const MultiIndex&M, int i)
00520         {
00521                 if (i==0) return M;
00522 
00523                 if (i<0)
00524                         return M + MultiIndex::Axis(-i+1); 
00525                 else
00526                         return M - MultiIndex::Axis( i-1);
00527         }
00528 
00529 protected:
00530 
00531 template <class Functor, int SuperDims>
00532 static  bool ForEachRecursion(Functor&F, const MultiIndex<SuperDims>&SuperIndex,
00533                               const MultiIndex&Start, const MultiIndex&End, const MultiIndex&Increment)
00534         {
00535                 for(index_t idx = Start.idx; idx < End.idx; idx += Increment.idx)
00536                 {
00537                 MultiIndex<SuperDims+1> SupraIndex(SuperIndex, idx); 
00538                         if (!MultiIndex<Dims-1>::ForEachRecursion(F,SupraIndex, Start, End, Increment) )
00539                                 return false;
00540                 } 
00541                 return true;
00542         }
00543 
00544 public:
00545 
00546 template <class Functor>
00547 static  bool ForEach(Functor&F, const MultiIndex&Start, const MultiIndex&End, const MultiIndex&Increment = MultiIndex(1) )
00548         {
00549                 for(index_t idx = Start.idx; idx < End.idx; idx += Increment.idx)
00550                 {
00551                 MultiIndex<1> SupraIndex(idx); 
00552                         if (!MultiIndex<Dims-1>::ForEachRecursion(F, SupraIndex, Start, End, Increment) )
00553                                 return false;
00554                 } 
00555                 return true;
00556         }
00557 };
00558 
00562 template <>
00563         class MultiIndex<1>
00564 {
00565         index_t idx;
00566 
00567 protected:
00568         typedef index_t Index;
00569         typedef Index value_type;
00570 
00571         MultiIndex(const MultiIndex&M, const MultiIndex&D, const Add&) throw()
00572         : idx(M.idx + D.maxidx() )
00573         {}
00574 
00575         MultiIndex(const MultiIndex&M, const MultiIndex&D, const Sub&) throw()
00576         : idx(M.idx - D.maxidx() )
00577         {}
00578 
00579         MultiIndex(const MultiIndex&M, const MultiIndex&D, const Mult&) throw()
00580         : idx(M.idx * D.maxidx() )
00581         {}
00582 
00583         MultiIndex(const MultiIndex&M, const MultiIndex&D, const Div&) throw()
00584         : idx(M.idx / D.maxidx() )
00585         {}
00586 
00587         MultiIndex(unsigned int bits, const ::Eagle::BinaryAnd&) throw()
00588         : idx( (bits & 1) ? 1 : 0 )
00589         {}
00590 
00591         MultiIndex(const MultiIndex&M, const MultiIndex&D, const ::Eagle::BinaryAnd&) throw()
00592         : idx( M.idx > D.maxidx() ? M.idx - D.maxidx() : D.maxidx() - M.idx )
00593         {}
00594 
00595         MultiIndex(const MultiIndex&M, const Power2Alignment&) throw()
00596         : idx( RoundUpToNextHighestPowerOf2( M.idx) )
00597         {}
00598 
00599 public:
00601         MultiIndex() throw()
00602         : idx(0)
00603         {}
00604 
00606         MultiIndex(const Index&I) throw()
00607         : idx(I)
00608         {}
00609 
00611         MultiIndex(const MultiIndex&Midx, const Index&I) throw()
00612         : idx(I)
00613         {}
00614 
00618         MultiIndex(const Index&lin_idx, const MultiIndex&Dimens, const CreateFromLinear&)
00619         : idx(lin_idx % Dimens.maxidx() )
00620         {}
00621 
00622 static  int     log2(index_t N) throw()
00623         {
00624                 if (N&1) return  0;
00625                 else     return -1;
00626         }
00627 
00628         int getSignedOrientation(const MultiIndex&B) const
00629         {
00630                 if (idx<B.idx) return +1;
00631                 if (idx>B.idx) return -1;
00632                 return 0;
00633         }
00634 
00635         const Index&operator[](int) const throw()
00636         {       return idx;      }
00637 
00638         Index&operator[](int) throw()
00639         {       return idx;      }
00640 
00641 
00642         const Index&maxidx() const throw()
00643         {       return idx; }
00644 
00645         Index&maxidx() throw()
00646         {       return idx;             }
00647 
00648         Index size() const throw()
00649         {       return idx;             }
00650 
00651 
00652         MultiIndex&operator+=(const MultiIndex&D) throw()
00653         {
00654                 idx += D.idx;
00655                 return *this;
00656         }
00657 
00658         MultiIndex&operator-=(const MultiIndex&D) throw()
00659         {
00660                 idx -= D.idx;
00661                 return *this;
00662         }
00663 
00664         MultiIndex operator+(const MultiIndex&D) throw()
00665         {
00666                 return MultiIndex(idx + D.idx);
00667         }
00668 
00669 
00670 
00671         friend MultiIndex clamp(const MultiIndex&M, const MultiIndex&D) throw()
00672         {
00673         MultiIndex retval; 
00674                 retval.idx      = clamp(M.idx, D.idx);
00675                 return retval;
00676         }
00677 
00678         friend MultiIndex clamp_m1(const MultiIndex&M, const MultiIndex&D) throw()
00679         {
00680         MultiIndex retval; 
00681                 retval.idx      = clamp(M.idx, D.idx-1);
00682                 return retval;
00683         }
00684 
00685 
00686         bool operator==(const MultiIndex&D) const throw()
00687         {       return idx == D.idx;    }
00688 
00690         bool operator!=(const MultiIndex&D) const throw()
00691         {
00692                 return idx != D.maxidx();
00693         }
00694 
00695         bool    inc(const MultiIndex&Dimens, const MultiIndex&Increment) throw()
00696         {       
00697                 idx+= Increment.maxidx();
00698                 if (idx >= Dimens.maxidx() )
00699                 {
00700                         idx = Index(0); 
00701                         return false;
00702                 } 
00703                 return true;
00704         }
00705         
00706         bool    inc(const MultiIndex&Dimens) throw()
00707         {       
00708                 idx++;
00709                 if (idx >= Dimens.maxidx() )
00710                 {
00711                         idx = Index(0); 
00712                         return false;
00713                 } 
00714                 return true;
00715         }
00716 
00717         Index   linear(const MultiIndex&/*Dimens*/) const throw()
00718         {       return idx;                             }
00719 
00720         bool    isWithin(const MultiIndex&Range) const throw()
00721         {
00722 //              return maxidx()>=0 && maxidx() < Range.maxidx();
00723                 return                maxidx() < Range.maxidx();
00724         }
00725 
00726 //               MSVC                           GCC             G++                   SGI C++
00727 #if     defined(_IOSTREAM_) || defined(_GLIBCXX_OSTREAM) || _CPP_IOSTREAM || defined(__SGI_STL_IOSTREAM)
00728 
00729         std::ostream&print(std::ostream&os, char sep) const
00730         //{     return os << IntValueTrait<Index>::getValue(idx);       }
00731         {       return os << idx;       }
00732 
00734         friend std::ostream&operator<<(std::ostream&os, const MultiIndex&MI)
00735         {
00736                 return MI.print(os, 'x');
00737         }
00738 #endif
00739 
00740 
00741 template <class Functor, int SuperDims>
00742 static  bool ForEachRecursion(Functor&F, const MultiIndex<SuperDims>&SuperIndex,
00743                               const MultiIndex&Start, const MultiIndex&End, const MultiIndex&Increment)
00744         {
00745                 for(index_t idx = Start.idx; idx < End.idx; idx += Increment.idx)
00746                 {
00747                 MultiIndex<SuperDims+1> SupraIndex(SuperIndex, idx); 
00748                         if (!F(SupraIndex) )
00749                                 return false;
00750                 } 
00751                 return true;
00752         }
00753 };
00754 
00755 
00756 inline MultiIndex<1> MIndex(index_t i0)
00757 {
00758 MultiIndex<1> m = i0;
00759         return m;
00760 }
00761 
00762 inline MultiIndex<2> MIndex(index_t i0, index_t i1)
00763 {
00764 MultiIndex<2> m; m = i0, i1;
00765         return m;
00766 }
00767 
00768 inline MultiIndex<3> MIndex(index_t i0, index_t i1, index_t i2)
00769 {
00770 MultiIndex<3> m; m = i0, i1, i2;
00771         return m;
00772 }
00773 
00774 inline MultiIndex<4> MIndex(index_t i0, index_t i1, index_t i2, index_t i3)
00775 {
00776 MultiIndex<4> m; m = i0, i1, i2, i3;
00777         return m;
00778 }
00779 
00780 
00781 
00782 
00783 
00784 
00785 //               MSVC                           GCC             G++                   SGI C++
00786 #if     defined(_IOSTREAM_) || defined(_GLIBCXX_OSTREAM) || _CPP_IOSTREAM || defined(__SGI_STL_IOSTREAM)
00787 template <int Dims>
00788   inline std::ostream&operator<<(std::ostream&os, const MultiIndex<Dims>&M)
00789 {
00790         return M.print( os<<'[', ',') << ']' ;
00791 }
00792 #endif
00793 
00794 
00799 template <int Dims>
00800 inline FixedArray<double, Dims> div(const MultiIndex<Dims>&Divisor, const MultiIndex<Dims>&Dividend)
00801 {
00802 FixedArray<double, Dims> result; 
00803 
00804         for(int i=0; i<Dims; i++)
00805                 result[i] = double(Divisor[i])/Dividend[i]; 
00806 
00807         return result;
00808 }
00809 
00810 
00811 #if 0
00812 
00816 template <int Dims>
00817         inline FixedArray<double, Dims> Modulo(double*a, double range=1.0)
00818         {
00819         FixedArray<double, Dims> m;
00820                 for(int i=0; i<Dims; i++)
00821                         m[i] = fmod(a[i] , range);
00822                 return m;
00823         }
00824 
00825 
00830 template <int Dims>
00831         inline FixedArray<double, Dims> Modulo(double*a, double range=1.0)
00832         {
00833         FixedArray<double, Dims> m;
00834                 for(int i=0; i<Dims; i++)
00835                         m[i] = fmod(a[i] , range);
00836                 return m;
00837         }
00838 
00843 template <int Dims>
00844         inline FixedArray<double, Dims> Modulo(const FixedArray<double, Dims>&a, const FixedArray<double, Dims>&range) 
00845         {
00846         MultiIndex<Dims, double> m; 
00847                 for(int i=0; i<Dims; i++)
00848                         m[i] = fmod(a[i] , range[i]);
00849                 return m;
00850         }
00851 #endif
00852 
00853 
00854 } // namespace Fiber
00855 
00856 #endif  // __Fiber_MultiIndex_HPP