MultiIndex.hpp

#ifndef __Fiber_MultiIndex_HPP
#define __Fiber_MultiIndex_HPP

#include <eagle/vecmath.hpp>
#include <eagle/Assignment.hpp>
#include <eagle/Operator.hpp>
#include <eagle/FixedArray.hpp>
#include "Iterator.hpp"

namespace Fiber
{

struct  CreateFromLinear {};

struct  Power2Alignment {};

template <int Dims>
        class MultiIndex : public MultiIndex<Dims-1>
{
        typedef MultiIndex<Dims-1>  base_t;

        index_t idx;

public:
        typedef index_t value_type;

        MultiIndex(const MultiIndex&M, const MultiIndex&D, const Add&) throw()
        : base_t(M,D, Add() )
        , idx(M.idx + D.maxidx())
        {}

        MultiIndex(const MultiIndex&M, const MultiIndex&D, const Sub&) throw()
        : base_t(M,D, Sub() )
        , idx(M.idx - D.maxidx())
        {}

        MultiIndex(const MultiIndex&M, const MultiIndex&D, const Mult&) throw()
        : base_t(M,D, Mult() )
        , idx(M.idx * D.maxidx())
        {}

        MultiIndex(const MultiIndex&M, const MultiIndex&D, const Div&)
        : base_t(M,D, Div() )
        , idx(M.idx / D.maxidx())
        {}

        MultiIndex(unsigned int bits, const ::Eagle::BinaryAnd&) throw()
        : base_t(bits, ::Eagle::BinaryAnd() )
        , idx( bits & (1<<(Dims-1)) ? 1 : 0 )
        {}

        MultiIndex(const MultiIndex&M, const MultiIndex&D, const ::Eagle::BinaryAnd&) throw()
        : base_t(M,D, ::Eagle::BinaryAnd() )
        , idx( M.idx > D.maxidx() ? M.idx - D.maxidx() : D.maxidx() - M.idx )
        {}

        MultiIndex(const MultiIndex&M, const Power2Alignment&) throw()
        : base_t(M, Power2Alignment() )
        , idx( RoundUpToNextHighestPowerOf2( M.idx) )
        {}

public:
        friend MultiIndex distance(const MultiIndex&M0, const MultiIndex&M1)
        {
                return MultiIndex(M0,M1, ::Eagle::BinaryAnd() );
        }

        friend MultiIndex bitalign(const MultiIndex&M) throw()
        {
                return MultiIndex(M, Power2Alignment() );
        }

static  MultiIndex BitIndex(unsigned int bits) throw()
        {
                return MultiIndex(bits, ::Eagle::BinaryAnd() );
        }

static  MultiIndex Axis(unsigned int orientation) throw()
        {
                return BitIndex( 1 << orientation );
        }

        index_t BitIndex() const throw()
        {
                return linear( MultiIndex(2) );
        }

static  int     log2(index_t N) throw()
        {
                if (N & (1<<(Dims-1) ) )
                    return Dims-1;

                return base_t::log2(N);
        }

        index_t Orientation() const throw()
        {
                return log2( BitIndex() );
        }       

        enum
        {
                dims = Dims,

                SIZE = Dims
        };

        int getSignedOrientation(const MultiIndex&B) const
        {
                if (idx<B.idx) return +Dims;
                if (idx>B.idx) return -Dims;
                return subidx().getSignedOrientation(B.subidx() );
        }

        MultiIndex() throw()
        : idx(0)
        {}

        explicit MultiIndex(const index_t&I) throw()
        : base_t(I), idx(I)
        {}

        MultiIndex(const MultiIndex<Dims-1>&Midx, const index_t&I) throw()
        : base_t(Midx), idx(I)
        {}

        MultiIndex(const MultiIndex<Dims-1>&m, const index_t&Slice, int SliceDirection)
        {
                for(int i=0, k=0; i<Dims; i++)
                {
                        if (i==SliceDirection)
                                (*this)[i] = Slice; 
                        else
                        {
                                (*this)[i] = m[k]; 
                                k++;
                        }
                }
        }

        MultiIndex(const index_t&lin_idx, const MultiIndex&Dimens, const CreateFromLinear&) throw()
        : base_t( lin_idx, Dimens, CreateFromLinear() )
        , idx( (lin_idx / Dimens.subidx().size()) % Dimens.maxidx() )
        {}

        const MultiIndex<Dims-1>&subidx() const throw()
        {       return *this;   }

        MultiIndex<Dims-1>&subidx() throw()
        {       return *this;   }

        const index_t&operator[](int i) const
        {       return (i>=Dims-1) ? idx : subidx()[i];         }

        index_t&operator[](int i) throw()
        {       return (i>=Dims-1) ? idx : subidx()[i];         }


        Eagle::Assignment<MultiIndex,0> operator=(const index_t&i)
        {
                return Eagle::Assignment<MultiIndex, 0>(*this, i);
        }

        const index_t&maxidx() const throw()
        {       return idx; }

        index_t&maxidx() throw()
        {       return idx; }

        index_t size() const throw()
        {
                return idx*subidx().size();     
        }

        bool operator!=(const MultiIndex&D) const throw()
        {
                return idx != D.maxidx() || subidx() != D.subidx() ;
        }

        bool operator==(const MultiIndex&M) const throw()
        {
                return idx == M.maxidx() && subidx() == M.subidx() ;
        }

        MultiIndex&operator+=(const MultiIndex&D) throw()
        {
                idx += D.maxidx(); 
                subidx() += D.subidx(); 
                return *this;
        }

        MultiIndex&operator-=(const MultiIndex&D) throw()
        {
                idx -= D.maxidx(); 
                subidx() -= D.subidx();
                return *this;
        }

/*
        friend MultiIndex clamp(const MultiIndex&M, const MultiIndex&D)
        {
        MultiIndex retval; 
                retval.idx      = clamp(M.idx     , D.idx);
                retval.subidx() = clamp(M.subidx(), D.subidx()); 
                return retval;
        }

        friend MultiIndex clamp_m1(const MultiIndex&M, const MultiIndex&D)
        {
        MultiIndex retval;
                retval.idx      = clamp_m1(M.idx     , D.idx-1);
                retval.subidx() = clamp_m1(M.subidx(), D.subidx());
                return retval;
        }
*/

        bool    inc(const MultiIndex&Dimens) throw()
        {
                idx ++;
                if ( idx >= Dimens.maxidx() )
                {
                bool    inc = subidx().inc(Dimens);
                        if (inc)
                                idx = index_t(0); 
                        return inc;
                } 
                return true;
        }

        bool    inc(const MultiIndex&Dimens, const MultiIndex&Increment) throw()
        {
                idx += Increment.maxidx();
                if ( idx >= Dimens.maxidx() )
                {
                bool    inc = subidx().inc( Dimens, Increment );
                        if (inc)
                                idx = index_t(0); 
                        return inc;
                } 
                return true;
        }

        index_t linear(const MultiIndex&Dimens) const throw()
        {       
//              return idx + Dimens.maxidx() * subidx().linear( Dimens.subidx() ); 
//              return idx * Dimens.maxidx() + subidx().linear( Dimens.subidx() ); 
//              return idx*subidx().size(); 

                assert(idx >= 0 && idx < Dimens.maxidx() );

                /*
                 Pretty sure, this formula can be optimized.
                 The size() call on the subdimension recursively multiplies
                 its elements, but this could be moved into the summation/
                 multiplication recursive call.
                */
        index_t L = Dimens.subidx().size();
                return idx * L + subidx().linear( Dimens.subidx() ); 
        }


//               MSVC                           GCC             G++                   SGI C++
#if     defined(_IOSTREAM_) || defined(_GLIBCXX_OSTREAM) || _CPP_IOSTREAM || defined(__SGI_STL_IOSTREAM)

        std::ostream&print(std::ostream&os, char sep) const
        {       
                base_t::print(os, sep); 
                return os << sep << idx;
        }

        friend std::ostream&operator<<(std::ostream&os, const MultiIndex&MI)
        {
                return MI.print(os, 'x');
        }
#endif

        bool    isWithin(const MultiIndex&Range) const throw()
        {
        // return maxidx()>=0 && maxidx() < Range.maxidx() &&
                return  maxidx() < Range.maxidx() &&
                        subidx().isWithin( Range );
        }

        bool    operator<(const MultiIndex&Range) const throw()
        {
                return isWithin( Range );
        }

        bool    operator>(const MultiIndex&Range) const throw()
        {
                return !isWithin( Range );
        }

        MultiIndex operator+(const MultiIndex&D) const throw()
        {       return MultiIndex(*this, D, Add() );    }

        MultiIndex operator-(const MultiIndex&D) const throw()
        {       return MultiIndex(*this, D, Sub() );    }

        MultiIndex operator*(const MultiIndex&D) const throw()
        {       return MultiIndex(*this, D, Mult() );   }

        MultiIndex operator/(const MultiIndex&D) const
        {       return MultiIndex(*this, D, Div() );    }


        MultiIndex operator+(int i) const
        {
        const MultiIndex&M = *this;
                if (i==0) return M;

                if (i<0)
                        return M - MultiIndex::Axis(-i+1); 
                else
                        return M + MultiIndex::Axis( i-1);
        }

        friend MultiIndex operator-(const MultiIndex&M, int i)
        {
                if (i==0) return M;

                if (i<0)
                        return M + MultiIndex::Axis(-i+1); 
                else
                        return M - MultiIndex::Axis( i-1);
        }

protected:

template <class Functor, int SuperDims>
static  bool ForEachRecursion(Functor&F, const MultiIndex<SuperDims>&SuperIndex,
                              const MultiIndex&Start, const MultiIndex&End, const MultiIndex&Increment)
        {
                for(index_t idx = Start.idx; idx < End.idx; idx += Increment.idx)
                {
                MultiIndex<SuperDims+1> SupraIndex(SuperIndex, idx); 
                        if (!MultiIndex<Dims-1>::ForEachRecursion(F,SupraIndex, Start, End, Increment) )
                                return false;
                } 
                return true;
        }

public:

template <class Functor>
static  bool ForEach(Functor&F, const MultiIndex&Start, const MultiIndex&End, const MultiIndex&Increment = MultiIndex(1) )
        {
                for(index_t idx = Start.idx; idx < End.idx; idx += Increment.idx)
                {
                MultiIndex<1> SupraIndex(idx); 
                        if (!MultiIndex<Dims-1>::ForEachRecursion(F, SupraIndex, Start, End, Increment) )
                                return false;
                } 
                return true;
        }
};

template <>
        class MultiIndex<1>
{
        index_t idx;

protected:
        typedef index_t Index;
        typedef Index value_type;

        MultiIndex(const MultiIndex&M, const MultiIndex&D, const Add&) throw()
        : idx(M.idx + D.maxidx() )
        {}

        MultiIndex(const MultiIndex&M, const MultiIndex&D, const Sub&) throw()
        : idx(M.idx - D.maxidx() )
        {}

        MultiIndex(const MultiIndex&M, const MultiIndex&D, const Mult&) throw()
        : idx(M.idx * D.maxidx() )
        {}

        MultiIndex(const MultiIndex&M, const MultiIndex&D, const Div&) throw()
        : idx(M.idx / D.maxidx() )
        {}

        MultiIndex(unsigned int bits, const ::Eagle::BinaryAnd&) throw()
        : idx( (bits & 1) ? 1 : 0 )
        {}

        MultiIndex(const MultiIndex&M, const MultiIndex&D, const ::Eagle::BinaryAnd&) throw()
        : idx( M.idx > D.maxidx() ? M.idx - D.maxidx() : D.maxidx() - M.idx )
        {}

        MultiIndex(const MultiIndex&M, const Power2Alignment&) throw()
        : idx( RoundUpToNextHighestPowerOf2( M.idx) )
        {}

public:
        MultiIndex() throw()
        : idx(0)
        {}

        MultiIndex(const Index&I) throw()
        : idx(I)
        {}

        MultiIndex(const MultiIndex&Midx, const Index&I) throw()
        : idx(I)
        {}

        MultiIndex(const Index&lin_idx, const MultiIndex&Dimens, const CreateFromLinear&)
        : idx(lin_idx % Dimens.maxidx() )
        {}

static  int     log2(index_t N) throw()
        {
                if (N&1) return  0;
                else     return -1;
        }

        int getSignedOrientation(const MultiIndex&B) const
        {
                if (idx<B.idx) return +1;
                if (idx>B.idx) return -1;
                return 0;
        }

        const Index&operator[](int) const throw()
        {       return idx;      }

        Index&operator[](int) throw()
        {       return idx;      }


        const Index&maxidx() const throw()
        {       return idx; }

        Index&maxidx() throw()
        {       return idx;             }

        Index size() const throw()
        {       return idx;             }


        MultiIndex&operator+=(const MultiIndex&D) throw()
        {
                idx += D.idx;
                return *this;
        }

        MultiIndex&operator-=(const MultiIndex&D) throw()
        {
                idx -= D.idx;
                return *this;
        }

        MultiIndex operator+(const MultiIndex&D) throw()
        {
                return MultiIndex(idx + D.idx);
        }



        friend MultiIndex clamp(const MultiIndex&M, const MultiIndex&D) throw()
        {
        MultiIndex retval; 
                retval.idx      = clamp(M.idx, D.idx);
                return retval;
        }

        friend MultiIndex clamp_m1(const MultiIndex&M, const MultiIndex&D) throw()
        {
        MultiIndex retval; 
                retval.idx      = clamp(M.idx, D.idx-1);
                return retval;
        }


        bool operator==(const MultiIndex&D) const throw()
        {       return idx == D.idx;    }

        bool operator!=(const MultiIndex&D) const throw()
        {
                return idx != D.maxidx();
        }

        bool    inc(const MultiIndex&Dimens, const MultiIndex&Increment) throw()
        {       
                idx+= Increment.maxidx();
                if (idx >= Dimens.maxidx() )
                {
                        idx = Index(0); 
                        return false;
                } 
                return true;
        }
        
        bool    inc(const MultiIndex&Dimens) throw()
        {       
                idx++;
                if (idx >= Dimens.maxidx() )
                {
                        idx = Index(0); 
                        return false;
                } 
                return true;
        }

        Index   linear(const MultiIndex&/*Dimens*/) const throw()
        {       return idx;                             }

        bool    isWithin(const MultiIndex&Range) const throw()
        {
//              return maxidx()>=0 && maxidx() < Range.maxidx();
                return                maxidx() < Range.maxidx();
        }

//               MSVC                           GCC             G++                   SGI C++
#if     defined(_IOSTREAM_) || defined(_GLIBCXX_OSTREAM) || _CPP_IOSTREAM || defined(__SGI_STL_IOSTREAM)

        std::ostream&print(std::ostream&os, char sep) const
        //{     return os << IntValueTrait<Index>::getValue(idx);       }
        {       return os << idx;       }

        friend std::ostream&operator<<(std::ostream&os, const MultiIndex&MI)
        {
                return MI.print(os, 'x');
        }
#endif


template <class Functor, int SuperDims>
static  bool ForEachRecursion(Functor&F, const MultiIndex<SuperDims>&SuperIndex,
                              const MultiIndex&Start, const MultiIndex&End, const MultiIndex&Increment)
        {
                for(index_t idx = Start.idx; idx < End.idx; idx += Increment.idx)
                {
                MultiIndex<SuperDims+1> SupraIndex(SuperIndex, idx); 
                        if (!F(SupraIndex) )
                                return false;
                } 
                return true;
        }
};


inline MultiIndex<1> MIndex(index_t i0)
{
MultiIndex<1> m = i0;
        return m;
}

inline MultiIndex<2> MIndex(index_t i0, index_t i1)
{
MultiIndex<2> m; m = i0, i1;
        return m;
}

inline MultiIndex<3> MIndex(index_t i0, index_t i1, index_t i2)
{
MultiIndex<3> m; m = i0, i1, i2;
        return m;
}

inline MultiIndex<4> MIndex(index_t i0, index_t i1, index_t i2, index_t i3)
{
MultiIndex<4> m; m = i0, i1, i2, i3;
        return m;
}






//               MSVC                           GCC             G++                   SGI C++
#if     defined(_IOSTREAM_) || defined(_GLIBCXX_OSTREAM) || _CPP_IOSTREAM || defined(__SGI_STL_IOSTREAM)
template <int Dims>
  inline std::ostream&operator<<(std::ostream&os, const MultiIndex<Dims>&M)
{
        return M.print( os<<'[', ',') << ']' ;
}
#endif


template <int Dims>
inline FixedArray<double, Dims> div(const MultiIndex<Dims>&Divisor, const MultiIndex<Dims>&Dividend)
{
FixedArray<double, Dims> result; 

        for(int i=0; i<Dims; i++)
                result[i] = double(Divisor[i])/Dividend[i]; 

        return result;
}


#if 0

template <int Dims>
        inline FixedArray<double, Dims> Modulo(double*a, double range=1.0)
        {
        FixedArray<double, Dims> m;
                for(int i=0; i<Dims; i++)
                        m[i] = fmod(a[i] , range);
                return m;
        }


template <int Dims>
        inline FixedArray<double, Dims> Modulo(double*a, double range=1.0)
        {
        FixedArray<double, Dims> m;
                for(int i=0; i<Dims; i++)
                        m[i] = fmod(a[i] , range);
                return m;
        }

template <int Dims>
        inline FixedArray<double, Dims> Modulo(const FixedArray<double, Dims>&a, const FixedArray<double, Dims>&range) 
        {
        MultiIndex<Dims, double> m; 
                for(int i=0; i<Dims; i++)
                        m[i] = fmod(a[i] , range[i]);
                return m;
        }
#endif


} // namespace Fiber

#endif  // __Fiber_MultiIndex_HPP