Shark Machine Learning Library

• About Shark
• Sourceforge
  • Project Summary
  • Downloads
  • Subversion Repository
• Getting Started
• Tutorials
• FAQ
• Main Modules
  • ReClaM
  • EALib
  • MOO-EALib
  • Fuzzy
• Tools
  • Mixture
  • Array
  • Rng
  • LinAlg
  • FileUtil

Array.h

Go to the documentation of this file.

//===========================================================================
//===========================================================================

#define __ARRAY_NO_GENERIC_IOSTREAM

#ifndef __ARRAY_H
#define __ARRAY_H


// backwards compatibility
#define Array2D Array
#define Array2DReference ArrayReference


#include <SharkDefs.h>
#include <algorithm>
#include <vector>
#include <typeinfo>


//==========================================================================
//
// doxygen declarations begin
//
//==========================================================================

//==========================================================================
//
// doxygen declarations end
//
//==========================================================================


// forward declaration
template< class T >
class Array;

// forward declaration
class ArrayIndex;

// forward declaration
template< class T >
class ArrayReference;

//===========================================================================
class ArrayBase
{
public:  // interface of public members

    //========================================================================
    unsigned ndim() const
    {
        return nd;
    }


    //========================================================================
    unsigned nelem() const
    {
        return ne;
    }


    //========================================================================
    bool samedim(const ArrayBase& v) const
    {
        if (nd == v.nd && ne == v.ne)
        {
            unsigned i;
            for (i = 0; i < nd && d[ i ] == v.d[ i ]; i++);
            return i == nd;
        }
        return false;
    }

    //========================================================================
    void resize(unsigned i, bool copy = false)
    {
        unsigned _d[ 1 ];
        _d[ 0 ] = i;
        resize_i(_d, 1, copy);
    }

    //========================================================================
    void resize(unsigned i, unsigned j, bool copy = false)
    {
        unsigned _d[ 2 ];
        _d[ 0 ] = i;
        _d[ 1 ] = j;
        resize_i(_d, 2, copy);
    }

    //========================================================================
    void resize(unsigned i, unsigned j, unsigned k, bool copy = false)
    {
        unsigned _d[ 3 ];
        _d[ 0 ] = i;
        _d[ 1 ] = j;
        _d[ 2 ] = k;
        resize_i(_d, 3, copy);
    }

    //========================================================================
    void resize(const std::vector< unsigned >& i, bool copy = false)
    {
        unsigned* _d = new unsigned[ i.size()];

        for (unsigned j = 0; j < i.size(); ++j)
            _d[ j ] = i[ j ];

        resize_i(_d, i.size(), copy);

        delete[ ] _d;
    }

    //========================================================================
    void resize(const ArrayBase& v, bool copy = false)
    {
        resize_i(v.d, v.nd, copy);
    }

    //========================================================================
    unsigned dim(unsigned i) const
    {
        RANGE_CHECK(i < nd)
        return d[ i ];
    }

    inline unsigned int rows() const
    {
        if (ndim() == 0) return 0;
        else return dim(0);
    }

    inline unsigned int cols() const
    {
        if (ndim() == 0) return 0;
        else return dim(ndim() - 1);
    }

    //========================================================================
    unsigned* dimvec()
    {
        return d;
    }

    //========================================================================
    const unsigned* dimvec() const
    {
        return d;
    }

    //========================================================================
    virtual ArrayBase* clone() const = 0;

    //========================================================================
    virtual ArrayBase* empty() const = 0;

    //========================================================================
    virtual ~ArrayBase()
    {
        if (! stat && nd)
        {
            delete[ ] d;
        }
    }

protected:  // interface of protected members (overloaded by array< T >)

    //========================================================================
    ArrayBase()
    {
        d    = 0;
        nd   = 0;
        ne   = 0;
        stat = false;
    }

    //=======================================================================
    ArrayBase(const ArrayBase& v)
    {
        nd    = v.nd;
        ne    = v.ne;
        stat  = false;

        if (nd)
        {
            d = new unsigned[ nd ];
            for (unsigned i = nd; i--; d[ i ] = v.d[ i ]);
        }
    }

    //=======================================================================
    virtual void resize_i(unsigned* d, unsigned nd, bool copy) = 0;

#ifndef __ARRAY_NO_GENERIC_IOSTREAM

    //=======================================================================
    virtual void readFrom(std::istream& is) = 0;

    //=======================================================================
    virtual void writeTo(std::ostream& os) const = 0;
#endif // !__ARRAY_NO_GENERIC_IOSTREAM

    //=======================================================================
    unsigned* d;

    unsigned  nd;

    unsigned  ne;

    bool      stat;

#ifndef __ARRAY_NO_GENERIC_IOSTREAM

    //=======================================================================
    friend inline std::istream& operator >> (std::istream& is, ArrayBase& a)
    {
        a.readFrom(is);
        IO_CHECK(is)
        return is;
    }

    //=======================================================================
    friend inline std::ostream& operator << (std::ostream& os, const ArrayBase& a)
    {
        a.writeTo(os);
        IO_CHECK(os)
        return os;
    }
#endif // __ARRAY_NO_GENERIC_IOSTREAM
};

//===========================================================================
template< class T >
class Array : public ArrayBase
{
    friend class ArrayReference< T >;
    friend class ArrayIndex;

public:

    typedef T* iterator;

    typedef T* const_iterator;


    //=======================================================================
    Array()
    {
        ne = 0;
        e = eEnd = 0;
    }


    //=======================================================================
    explicit Array(unsigned i)
    {
        d      = new unsigned[ 1 ];
        d[ 0 ] = i;
        nd     = 1;
        ne     = i;
        e      = 0;
        if (ne)
        {
            e = new T[ ne ];
        }
        stat = false;
        eEnd = e + ne;
    }

    //=======================================================================
    Array(unsigned i, unsigned j)
    {
        d      = new unsigned[ 2 ];
        d[ 0 ] = i;
        d[ 1 ] = j;
        nd     = 2;
        ne     = i * j;
        e      = 0;
        if (ne) e = new T[ ne ];
        stat = false;
        eEnd = e + ne;
    }

    //=======================================================================
    Array(unsigned i, unsigned j, unsigned k)
    {
        d      = new unsigned[ 3 ];
        d[ 0 ] = i;
        d[ 1 ] = j;
        d[ 2 ] = k;
        nd     = 3;
        ne     = i * j * k;
        e      = 0;
        if (ne)
        {
            e = new T[ ne ];
        }
        stat = false;
        eEnd = e + ne;
    }

    //=======================================================================
    Array(std::vector< T >& v)
    {
        d      = new unsigned[ 1 ];
        d[ 0 ] = v.size();
        nd     = 1;
        ne     = v.size();
        e      = 0;
        if (ne)
        {
            e = new T[ ne ];
            for (unsigned i = ne; i--; e[ i ] = v[ i ]);
        }
        stat = false;
        eEnd = e + ne;
    }

    //=======================================================================
    Array(Array< T >& v, bool)
    {
        d      = v.d;
        e      = v.e;
        nd     = v.nd;
        ne     = v.ne;
        stat   = v.stat;
        v.nd   = v.ne = 0;
        v.stat = true;
        eEnd   = e + ne;
    }

    //=======================================================================
    Array(const Array< T >& v) : ArrayBase(v)
    {
        if (ne)
        {
            e = new T[ ne ];
            for (unsigned i = ne; i--; e[ i ] = v.e[ i ]);
            eEnd = e + ne;
        }
        else
        {
            e = eEnd = 0;
        }
    }

    //=======================================================================
    ~Array()
    {
        if (! stat && ne)
        {
            if (ne)
            {
                delete[ ] e;
            }
        }
    }

    inline unsigned int rows() const
    {
        return ArrayBase::rows();
    }

    inline unsigned int cols() const
    {
        return ArrayBase::cols();
    }

    //=======================================================================
    T& operator()()
    {
        SIZE_CHECK(nd == 0)
        RANGE_CHECK(ne == 1)
        return e[ 0 ];
    }

    //=======================================================================
    const T& operator()() const
    {
        SIZE_CHECK(nd == 0)
        RANGE_CHECK(ne == 1)
        return e[ 0 ];
    }

    //=======================================================================
    T& operator()(unsigned i)
    {
        SIZE_CHECK(nd == 1)
        RANGE_CHECK(i < d[ 0 ])
        return e[ i ];
    }

    //=======================================================================
    const T& operator()(unsigned i) const
    {
        SIZE_CHECK(nd == 1)
        RANGE_CHECK(i < d[ 0 ])
        return e[ i ];
    }

    //=======================================================================
    T& operator()(unsigned i, unsigned j)
    {
        SIZE_CHECK(nd == 2)
        RANGE_CHECK(i < d[ 0 ])
        RANGE_CHECK(j < d[ 1 ])
        return e[ d[ 1 ] * i + j ];
    }

    //=======================================================================
    const T& operator()(unsigned i, unsigned j) const
    {
        SIZE_CHECK(nd == 2)
        RANGE_CHECK(i < d[ 0 ])
        RANGE_CHECK(j < d[ 1 ])
        return e[ d[ 1 ] * i + j ];
    }

    //=======================================================================
    T& operator()(unsigned i, unsigned j, unsigned k)
    {
        SIZE_CHECK(nd == 3)
        RANGE_CHECK(i < d[ 0 ])
        RANGE_CHECK(j < d[ 1 ])
        RANGE_CHECK(k < d[ 2 ])
        return e[ d[ 2 ] * d[ 1 ] * i + d[ 2 ] * j + k ];
    }

    //=======================================================================
    const T& operator()(unsigned i, unsigned j, unsigned k) const
    {
        SIZE_CHECK(nd == 3)
        RANGE_CHECK(i < d[ 0 ])
        RANGE_CHECK(j < d[ 1 ])
        RANGE_CHECK(k < d[ 2 ])
        return e[ d[ 2 ] * d[ 1 ] * i + d[ 2 ] * j + k ];
    }


    //=======================================================================
    T& operator()(const std::vector< unsigned >& i)
    {
        unsigned l, m, n;
        SIZE_CHECK(i.size() == nd)
#ifndef NDEBUG
        for (l = nd; l--;)
        {
            RANGE_CHECK(i[ l ] < d[ l ])
        }
#endif
        for (l = nd, m = 0, n = 1; l--;)
        {
            m += i[ l ] * n;
            n *= d[ l ];
        }

        return e[ m ];
    }


    //=======================================================================
    const T& operator()(const std::vector< unsigned >& i) const
    {
        unsigned l, m, n;
        SIZE_CHECK(i.size() == nd)
#ifndef NDEBUG
        for (l = nd; l--;)
        {
            RANGE_CHECK(i[ l ] < d[ l ])
        }
#endif
        for (l = nd, m = 0, n = 1; l--;)
        {
            m += i[ l ] * n;
            n *= d[ l ];
        }

        return e[ m ];
    }


    // Forward reference (see end of class ArrayReference for definition)
    ArrayReference< T > operator [ ](unsigned i);

    // Forward reference (see end of class ArrayReference for definition)
    const ArrayReference< T > operator [ ](unsigned i) const;

    //========================================================================
    Array< T >& operator = (const T& v)
    {
        for (unsigned i = ne; i--; e[ i ] = v);
        return *this;
    }




    //========================================================================
    Array< T >& operator = (const std::vector< T >& v)
    {
        resize(v.size());
        for (unsigned i = ne; i--; e[ i ] = v[ i ]);
        return *this;
    }



    //========================================================================
    Array< T >& operator = (const Array< T >& v)
    {
        //
        // handle the special case of assigning an empty array
        // in order to avoid uninitialized memory read (purify)
        //
        if (v.nd)
        {
            resize_i(v.d, v.nd, false);
            for (unsigned i = ne; i--; e[ i ] = v.e[ i ]);
        }
        else
        {
            if (nd)
            {
                delete[ ] d;
            }
            if (ne)
            {
                delete[ ] e;
            }
            nd = ne = 0;
        }

        return *this;
    }



    //========================================================================
    T& elem(unsigned i)
    {
        RANGE_CHECK(i < ne)
        return e[ i ];
    }


    //========================================================================
    const T& elem(unsigned i) const
    {
        RANGE_CHECK(i < ne)
        return e[ i ];
    }


    //========================================================================
    T* elemvec()
    {
        return e;
    }


    //========================================================================
    const T* elemvec() const
    {
        return e;
    }


    //========================================================================
    Array< T >& append_elem(const T& w)
    {
        SIZE_CHECK(nd == 0 || nd == 1)

        resize(ne + 1, true);
        e[ ne - 1 ] = w;
        return *this;
    }


    //========================================================================
    Array< T >& append_elems(const Array< T >& w)
    {
        SIZE_CHECK((nd == 0 || nd == 1) && w.nd == 1)

        unsigned i = 0, j = ne;

        resize(ne + w.ne, true);

        while (i < w.ne)
        {
            e[ j++ ] = w.e[ i++ ];
        }

        return *this;
    }



    //========================================================================
    Array< T >& append_rows(const Array< T >& y)
    {
        SIZE_CHECK(ndim() == 0         ||
                   ndim() == y.ndim() ||
                   ndim() == y.ndim() + 1)

        unsigned i, pos = nelem();
        Array< unsigned > da;//is used as to store the new dimensions

        // Append rows to empty array:
        if (ndim() == 0)
        {
            da.resize(y.ndim() + 1);
            da(0) = 1;
            for (i = 0; i < y.ndim(); ++i)
            {
                da(i + 1) = y.dim(i);
            }
        }
        // Append rows if dimensions are equal
        else if (ndim() == y.ndim())
        {
          //dimarr() returns a copy of the current arrays 
          //dimension vector as an array
            da = dimarr();
            da(0) = dim(0) + y.dim(0);
        }
        // Append rows if y's dimension is one less the 
        // the arrays dimension
        else if (ndim() == y.ndim() + 1)
        {
            da = dimarr();
            da(0) ++;
        }
        // Otherwise do nothing and return the original array
        else
        {
            return *this;
        }
        // resize the array, keep its elements and fill 
        // the new rows with zeros
        resize_i(da.elemvec(), da.nelem(), true);
        // add the new elements
        for (i = 0; i < y.nelem();)
        {
            elem(pos++) = y.elem(i++);
        }

        return *this;
    }




    //========================================================================
    Array< T >& remove_row(unsigned i)
    {
        RANGE_CHECK(ndim() > 0 && i < dim(0))

        unsigned k;
        unsigned size = nelem() / dim(0);// number of rows

        Array< T        > y;// will be the new array
        Array< unsigned > da(dimarr());// stores the new dimensions

        da(0)--;

        y.resize_i(da.elemvec(), da.nelem(), false);
        // copy all rows with indices smaller than the deleted row
        for (k = 0; k < i * size; ++k)
        {
            y.elem(k) = elem(k);
        }
        // copy all rows with indices greater than the deleted row
        for (; k < y.nelem(); ++k)
        {
            y.elem(k) = elem(k + size);
        }
        // replace the current array by y
        return *this = y;
    }




    //========================================================================
    Array< T > remove_col(unsigned k) const
    {
        RANGE_CHECK(ndim() > 0 && k < dim(ndim() - 1))

        unsigned i, j, xi, zi;
        unsigned xlast = dim(ndim() - 1);
        //number of columns
        unsigned num   = nelem() / xlast;
        //the new array
        Array< T > z;
        //stores the dimensions of the new array
        Array< unsigned > da(dimarr());
        //new array has one column less
        da(ndim() - 1) = xlast - 1;
        //resize z according to da
        z.resize_i(da.elemvec(), da.nelem(), false);
        // copy all kept columns in z 
        for (xi = zi = i = 0; i < num; ++i)
        {
          // copy all columns with a lower index than k
            for (j = 0; j < k; ++j)
            {
                z.elem(zi++) = elem(xi++);
            }
          // copy all columns with a greater index than k
            for (++xi, ++j; j < xlast; ++j)
            {
                z.elem(zi++) = elem(xi++);
            }
        }
        //return new array
        return Array< T > (z, true);
    }




    //========================================================================
    Array<T>& remove_cols(const Array<unsigned> idx)
    {
        RANGE_CHECK(ndim() > 0 && idx.ndim() == 1);

        if (idx.nelem() == 0) return *this;
        std::sort(idx.begin(), idx.end());

        T* source = e;
        T* dest = e;
        int i;
        int lastDim = ndim() - 1;
        int r, rc = idx.nelem();
        int f, fc = dim(lastDim);
        RANGE_CHECK(idx(rc - 1) < fc);
        for (i=0, r=0, f=0; i<ne; i++)
        {
            if (r < rc && f == idx(r))
            {
                source++;
                r++;
            }
            else
            {
                *dest = *source;
                source++;
                dest++;
            }
            f++;
            if (f == fc)
            {
                f = 0;
                r = 0;
            }
        }

        std::vector<unsigned> newdim(ndim());
        for (i=0; i<lastDim; i++) newdim[i] = d[i];
        newdim[lastDim] = fc - rc;
        resize_i(&newdim[0], ndim(), true);

        return *this;
    }




    //========================================================================
    Array< T > & remove_cols(unsigned i, unsigned j)
    {
        Array< unsigned > idx(2);

        idx(0) = i;
        idx(1) = j;

        return remove_cols(idx);
    }



    //========================================================================
    Array< T > & remove_cols(unsigned i, unsigned j, unsigned k)
    {
        Array< unsigned > idx(3);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;

        return remove_cols(idx);
    }



    //========================================================================
    Array< T > & remove_cols(unsigned i, unsigned j, unsigned k,
                           unsigned l)
    {
        Array< unsigned > idx(4);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;
        idx(3) = l;

        return remove_cols(idx);
    }



    //========================================================================
    Array< T > & remove_cols(unsigned i, unsigned j, unsigned k,
                           unsigned l, unsigned m)
    {
        Array< unsigned > idx(5);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;
        idx(3) = l;
        idx(4) = m;

        return remove_cols(idx);
    }



    //========================================================================
    Array< T > & remove_cols(unsigned i, unsigned j, unsigned k,
                           unsigned l, unsigned m, unsigned n)
    {
        Array< unsigned > idx(6);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;
        idx(3) = l;
        idx(4) = m;
        idx(5) = n;

        return remove_cols(idx);
    }



    //========================================================================
    Array< T > & remove_cols(unsigned i, unsigned j, unsigned k,
                           unsigned l, unsigned m, unsigned n,
                           unsigned o)
    {
        Array< unsigned > idx(7);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;
        idx(3) = l;
        idx(4) = m;
        idx(5) = n;
        idx(6) = o;

        return remove_cols(idx);
    }



    //========================================================================
    Array< T > & remove_cols(unsigned i, unsigned j, unsigned k,
                           unsigned l, unsigned m, unsigned n,
                           unsigned o, unsigned p)
    {
        Array< unsigned > idx(8);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;
        idx(3) = l;
        idx(4) = m;
        idx(5) = n;
        idx(6) = o;
        idx(7) = p;

        return remove_cols(idx);
    }


    //========================================================================
    Array< T > append_cols(const Array< T >& y) const
    {
        SIZE_CHECK(ndim() == 0 || ndim() == y.ndim())

        if (ndim() == 0)
        {
            return y;
        }
        else
        {
            unsigned i, j, xi, yi, zi;
            //number of last column
            unsigned xlast = dim(ndim() - 1);
            //number of last columns in y
            unsigned ylast = y.dim(y.ndim() - 1);
            unsigned num;
            //new array
            Array< T > z;
            //stores new dimensions
            Array< unsigned > da(dimarr());
            //new number of columns
            da(ndim() - 1) = xlast + ylast;
            //resize new array according to da
            z.resize_i(da.elemvec(), da.nelem(), false);

            num = z.nelem() / (xlast + ylast);//number of rows
            for (i = xi = yi = zi = 0; i < num; ++i)
            {
              //copy all elements from old array
                for (j = 0; j < xlast; ++j)
                {
                    z.elem(zi++) = elem(xi++);
                }
              //copy all elements from y
                for (j = 0; j < ylast; ++j)
                {
                    z.elem(zi++) = y.elem(yi++);
                }
            }
            //return new array
            return Array< T > (z, true);
        }
    }




    //========================================================================
    Array< T > subarr(unsigned from, unsigned to) const
    {
        SIZE_CHECK(ndim() > 0)
        RANGE_CHECK(from <= to && from < dim(0) && to < dim(0))

        unsigned i, j;
        //number of first element
        unsigned fromi = from   * nelem() / dim(0);
        //number of last element
        unsigned toi   = (to + 1) * nelem() / dim(0);
        //subarray
        Array< T > z;
        //stores dimensions of subarray
        Array< unsigned > da(dimarr());
        // number of rows in subarray
        da(0) = to - from + 1;
        //resize subarray
        z.resize_i(da.elemvec(), da.nelem(), false);
        //copy elements
        for (i = fromi, j = 0; i < toi; z.elem(j++) = elem(i++));
        //return subarray
        return Array< T > (z, true);
    }




    //========================================================================
    Array< unsigned > pos2idx(unsigned p)
    {
        RANGE_CHECK(p < ne)

        Array< unsigned > idx(ndim());
        for (unsigned i = 0, j = nelem(); i < ndim(); ++i)
        {
            p %= j;
            j /= dim(i);
            idx(i) = p / j;
        }
        return idx;
    }



    //========================================================================
    void pos2idx(unsigned p, Array<unsigned>& idx)
    {
        RANGE_CHECK(p < ne)

        idx.resize(ndim());
        for (unsigned i = 0, j = nelem(); i < ndim(); ++i)
        {
            p %= j;
            j /= dim(i);
            idx(i) = p / j;
        }
    }


    //========================================================================
    Array< T > row(unsigned i) const
    {
        return (*this)[ i ];
    }


    //========================================================================
    Array< T > rows(const Array< unsigned >& idx) const
    {
        SIZE_CHECK(ndim() > 0)
        SIZE_CHECK(idx.ndim() <= 1)

        unsigned i, j, k, l;
        //number of rows
        unsigned n = nelem() / dim(0);
        //subarray
        Array< T > z;
        //stores dimensions of subarray
        Array< unsigned > da(dimarr());

        da(0) = idx.nelem();
        //resize subarray
        z.resize_i(da.elemvec(), da.nelem(), false);
        //copy all elements mentioned in idx 
        for (k = i = 0; i < idx.nelem(); ++i)
        {
            for (l = idx(i) * n, j = 0; j < n; ++j)
            {
                z.elem(k++) = elem(l++);
            }
        }
        //return subarray
        return Array< T > (z, true);
    }



    //========================================================================
    Array< T > rows(unsigned i) const
    {
        Array< unsigned > idx(1);

        idx(0) = i;

        return rows(idx);
    }



    //========================================================================
    Array< T > rows(unsigned i, unsigned j) const
    {
        Array< unsigned > idx(2);

        idx(0) = i;
        idx(1) = j;

        return rows(idx);
    }



    //========================================================================
    Array< T > rows(unsigned i, unsigned j, unsigned k) const
    {
        Array< unsigned > idx(3);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;

        return rows(idx);
    }



    //========================================================================
    Array< T > rows(unsigned i, unsigned j, unsigned k,
                    unsigned l) const
    {
        Array< unsigned > idx(4);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;
        idx(3) = l;

        return rows(idx);
    }



    //========================================================================
    Array< T > rows(unsigned i, unsigned j, unsigned k,
                    unsigned l, unsigned m) const
    {
        Array< unsigned > idx(5);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;
        idx(3) = l;
        idx(4) = m;

        return rows(idx);
    }



    //========================================================================
    Array< T > rows(unsigned i, unsigned j, unsigned k,
                    unsigned l, unsigned m, unsigned n) const
    {
        Array< unsigned > idx(6);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;
        idx(3) = l;
        idx(4) = m;
        idx(5) = n;

        return rows(idx);
    }



    //========================================================================
    Array< T > rows(unsigned i, unsigned j, unsigned k,
                    unsigned l, unsigned m, unsigned n,
                    unsigned o) const
    {
        Array< unsigned > idx(7);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;
        idx(3) = l;
        idx(4) = m;
        idx(5) = n;
        idx(6) = o;

        return rows(idx);
    }



    //========================================================================
    Array< T > rows(unsigned i, unsigned j, unsigned k,
                    unsigned l, unsigned m, unsigned n,
                    unsigned o, unsigned p) const
    {
        Array< unsigned > idx(8);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;
        idx(3) = l;
        idx(4) = m;
        idx(5) = n;
        idx(6) = o;
        idx(7) = p;

        return rows(idx);
    }




    //========================================================================
    Array< T > col(unsigned i) const
    {
        SIZE_CHECK(ndim() > 0)
        RANGE_CHECK(i < dim(ndim() - 1))
        //subarray
        Array< T > z;
        //stores dimensions of subarray
        Array< unsigned > da(dimarr());
        //resize subarray
        z.resize_i(da.elemvec(), da.nelem() - 1, false);
        //copy i'th  col
        for (unsigned k = 0; k < z.nelem(); ++k, i += dim(ndim() - 1))
        {
            z.elem(k) = elem(i);
        }
        //return subarry
        return z;
    }



    //========================================================================
    Array< T > cols(const Array< unsigned >& idx) const
    {
        SIZE_CHECK(ndim() > 0)
        SIZE_CHECK(idx.ndim() == 1)

        unsigned j, k, m;
        //number of columns
        unsigned l = dim(ndim() - 1);
        //number of rows
        unsigned n = nelem() / l;
        //subarray
        Array< T > z;
        Array< unsigned > ix(idx);
        //stores dimensions of subarray
        Array< unsigned > da(dimarr());

        da(ndim() - 1) = ix.nelem();
        //resize subarray
        z.resize_i(da.elemvec(), da.nelem(), false);
        //copy all elements mentioned in idx
        for (k = m = 0; k < n; ++k)
        {
            for (j = 0; j < ix.nelem(); ++j, ++m)
            {
                z.elem(m) = elem(ix(j));
                ix(j) += l;//switch to next row
            }
        }
        //return subarray
        return z;
    }



    //========================================================================
    Array< T > cols(unsigned i) const
    {
        Array< unsigned > idx(1);

        idx(0) = i;

        return cols(idx);
    }



    //========================================================================
    Array< T > cols(unsigned i, unsigned j) const
    {
        Array< unsigned > idx(2);

        idx(0) = i;
        idx(1) = j;

        return cols(idx);
    }



    //========================================================================
    Array< T > cols(unsigned i, unsigned j, unsigned k) const
    {
        Array< unsigned > idx(3);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;

        return cols(idx);
    }



    //========================================================================
    Array< T > cols(unsigned i, unsigned j, unsigned k,
                    unsigned l) const
    {
        Array< unsigned > idx(4);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;
        idx(3) = l;

        return cols(idx);
    }



    //========================================================================
    Array< T > cols(unsigned i, unsigned j, unsigned k,
                    unsigned l, unsigned m) const
    {
        Array< unsigned > idx(5);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;
        idx(3) = l;
        idx(4) = m;

        return cols(idx);
    }



    //========================================================================
    Array< T > cols(unsigned i, unsigned j, unsigned k,
                    unsigned l, unsigned m, unsigned n) const
    {
        Array< unsigned > idx(6);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;
        idx(3) = l;
        idx(4) = m;
        idx(5) = n;

        return cols(idx);
    }



    //========================================================================
    Array< T > cols(unsigned i, unsigned j, unsigned k,
                    unsigned l, unsigned m, unsigned n,
                    unsigned o) const
    {
        Array< unsigned > idx(7);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;
        idx(3) = l;
        idx(4) = m;
        idx(5) = n;
        idx(6) = o;

        return cols(idx);
    }



    //========================================================================
    Array< T > cols(unsigned i, unsigned j, unsigned k,
                    unsigned l, unsigned m, unsigned n,
                    unsigned o, unsigned p) const
    {
        Array< unsigned > idx(8);

        idx(0) = i;
        idx(1) = j;
        idx(2) = k;
        idx(3) = l;
        idx(4) = m;
        idx(5) = n;
        idx(6) = o;
        idx(7) = p;

        return cols(idx);
    }



    //========================================================================
    Array< T >& transpose()
    {
        switch (ndim())
        {
        case 0 :
        case 1 :
            break;
        case 2 :
        {
            int d0 = d[1];
            int d1 = d[0];
            int i, j, k, s;
            for (i = 0, s = d0 * d1; s > 0; i++)
            {
                for (j = (i % d1) * d0 + i / d1; j > i; j = (j % d1) * d0 + j / d1);
                if (j < i) continue;
                for (k = i, j = (i % d1) * d0 + i / d1; j != i; k = j, j = (j % d1) * d0 + j / d1)
                {
                    std::swap(e[k], e[j]);
                    s--;
                }
                s--;
            }
            d[0] = d0;
            d[1] = d1;
            break;
        }
        default :
        {
            int d0 = d[1];
            int d1 = d[0];
            int i, j, k, s;
            unsigned a, b, c, size = nelem() / (d0 * d1);
            for (i = 0, s = d0 * d1; s > 0; i++)
            {
                for (j = (i % d1) * d0 + i / d1; j > i; j = (j % d1) * d0 + j / d1);
                if (j < i) continue;
                for (k = i, j = (i % d1) * d0 + i / d1; j != i; k = j, j = (j % d1) * d0 + j / d1)
                {
                    a = k * size;
                    b = j * size;
                    for (c = 0; c < size; c++)
                    {
                        std::swap(e[a], e[b]);
                        a++;
                        b++;
                    }
                    s--;
                }
                s--;
            }
            d[0] = d0;
            d[1] = d1;
            break;
        }
        }

        return *this;
    }

    //========================================================================
    Array<T>& rotate_rows(int n)
    {
        SIZE_CHECK(ndim() > 0);

        int f, fc = d[0];
        int i, ic = ne / fc;
        std::vector<T> tmp(fc);

        n = n % fc;
        if (n < 0) n += fc;
        if (n == 0) return *this;
        RANGE_CHECK(0 <= n && n < fc);

        T* data = e;
        for (i=0; i<ic; i++)
        {
            for (f=0; f<fc; f++) tmp[f] = data[ic * f];
            for (f=0; f<n; f++) data[ic * f] = tmp[f - n + fc];
            for (; f<fc; f++) data[ic * f] = tmp[f - n];
            data++;
        }

        return *this;
    }


    //========================================================================
    Array<T>& rotate_cols(int n)
    {
        SIZE_CHECK(ndim() > 0);

        int lastDim = nd - 1;
        int f, fc = d[lastDim];
        int i, ic = ne / fc;
        std::vector<T> tmp(fc);

        n = n % fc;
        if (n < 0) n += fc;
        if (n == 0) return *this;
        RANGE_CHECK(0 < n && n < fc);

        T* data = e;
        for (i=0; i<ic; i++)
        {
            for (f=0; f<fc; f++) tmp[f] = data[f];
            for (f=0; f<n; f++) data[f] = tmp[f - n + fc];
            for (; f<fc; f++) data[f] = tmp[f - n];
            data += fc;
        }

        return *this;
    }


    //========================================================================
    Array< unsigned > dimarr() const
    {
        Array< unsigned > dim(nd);
        for (unsigned i = nd; i--; dim(i) = d[ i ]);
        return dim;
    }



    //========================================================================
    ArrayBase* clone() const
    {
        return new Array< T > (*this);
    }


    //========================================================================
    ArrayBase* empty() const
    {
        return new Array< T >();
    }

#ifdef _WIN32
    bool operator == (const Array< T >&) const
    {
        return false;
    }

    bool operator != (const Array< T >&) const
    {
        return false;
    }

    bool operator < (const Array< T >&) const
    {
        return false;
    }

    bool operator > (const Array< T >&) const
    {
        return false;
    }
#endif


    //=======================================================================
    inline iterator begin()
    {
        return e;
    }


    //=======================================================================
    inline const const_iterator begin() const
    {
        return e;
    }

    //=======================================================================
    inline iterator end()
    {
        return eEnd;
    }

    //=======================================================================
    inline const const_iterator end()   const
    {
        return eEnd;
    }

protected:

    //=======================================================================
    T*        e;

    T*        eEnd;


    //=======================================================================
    void resize_i(unsigned* _d, unsigned _nd, bool copy)
    {
        unsigned j, _ne;

        // Adopt dimension vector if necessary:
        if (nd != _nd)
        {
            // Can only resize non-static arrays:
            SIZE_CHECK(! stat)

            // Delete dimension vector if present:
            if (nd)
            {
                delete[ ] d;
            }

            // Create new dimension vector:
            if ((nd = _nd) > 0)
            {
                d = new unsigned[ _nd ];
            }
        }

#ifndef NDEBUG
        if (stat)
        {
            for (j = _nd; j--;)
            {
                SIZE_CHECK(d[ j ] == _d[ j ])
            }
        }
#endif

        //calculate number of elements
        for (j = _nd, _ne = 1; j--; _ne *= (d[ j ] = _d[ j ]));
        if (_nd == 0)
        {
            _ne = 0;
        }

        if (ne != _ne)  //add or delete elements
        {
            SIZE_CHECK(! stat);
#ifndef NDEBUG
            if (nd != _nd) throw ("Array resize warning: copying elements will not arrange elements correctly");
            else
            {
                for (j = 1; j < _nd; j++)
                {
                    if (d[j] != _d[j])
                        throw ("Array resize warning: copying elements will not arrange elements correctly");
                }
            }
#endif

            if (copy)   //does not consider structure of elements
            {
                T* _e = e;  //create new pointer to old element vector
                if (_ne)    //create new element vector
                {
                    e = new T[ _ne ];
                }
                for (j = ne < _ne ? ne : _ne; j--; e[ j ] = _e[ j ]);//copy existing elements, j=min(ne,_ne)
                if (ne)//delete elements, array is empty now
                {
                    delete[ ] _e;
                }
            }
            else//without copying
            {
                if (ne)//delete element vector
                {
                    delete[ ] e;
                }
                if (_ne)//create new element vector
                {
                    e = new T[ _ne ];
                }
            }
            ne = _ne;//update number of elements in array
        }
        eEnd = e + ne;//update pointer to end of array
    }


    //=======================================================================
    Array(unsigned* _d, T* _e, unsigned _nd, unsigned _ne, bool _stat)
    {
        d    = _d;
        e    = _e;
        nd   = _nd;
        ne   = _ne;
        stat = _stat;
        eEnd = e + ne;
    }

#ifndef __ARRAY_NO_GENERIC_IOSTREAM


    //=======================================================================
    void readFrom(std::istream& is)
    {
        const unsigned     MaxLen = 1024;
        char               s[ MaxLen ];
        char*              p;
        unsigned           i;
        std::vector< unsigned > idx; // dimension vector

        // Skip leading whitespaces
        is >> std::ws;
        is.getline(s, MaxLen);

        // Get end of type definition:
        p = strchr(s, '>');
        // Go to begin of dimension sizes:
        if (p && *p) p++;
        if (p && *p) p++;

        // Get single dimension sizes, store them in vector
        // and resize current array:
        TYPE_CHECK(p && *p)
        if (p && *p)
        {
            do
            {
                if (*p == ',')
                {
                    p++;
                }
                idx.push_back(unsigned(strtol(p, &p, 10)));
            }
            while (p && *p == ',');

            TYPE_CHECK(p && *p == ')')

            resize(idx);

            // Get elements of the array:
            for (i = 0; i < ne && is.good(); i++)
            {
                is >> e[ i ];
            }
        }

        is >> std::ws;
    }



    //=======================================================================
    void writeTo(std::ostream& os) const
    {
        unsigned i;

        os << "Array<" << typeid(T).name() << ">(";
        for (i = 0; i < nd; i++)
        {
            if (i) os << ',';
            os << d[ i ];
        }
        os << ")\n";

        for (i = 0; i < ne; i++)
        {
            if (i) os << '\t';
            os << e[ i ];
        }
        os << std::endl;
    }
#endif // !__ARRAY_NO_GENERIC_IOSTREAM
};


//===========================================================================
template< class T >
class ArrayReference : public Array< T >
{
public:

    //========================================================================
    ArrayReference()
            : Array< T >()
    {
        this->stat = true;
        this->d = NULL;
        this->e = NULL;
        this->nd = 0;
        this->ne = 0;
    }


    //========================================================================
    ArrayReference(unsigned* _d, T* _e, unsigned _nd, unsigned _ne)
            : Array< T >()
    {
        this->d    = _d;
        this->e    = _e;
        this->nd   = _nd;
        this->ne   = _ne;
        this->stat = true;
        this->eEnd = this->e + this->ne;
    }


    //========================================================================
    ~ArrayReference()
    { }



    //========================================================================
    ArrayReference< T > operator = (const T& v)
    {
        Array< T >::operator = (v);
        return *this;
    }



    //========================================================================
    ArrayReference< T > operator = (const std::vector< T >& v)
    {
        Array< T >::operator = (v);
        return *this;
    }



    //========================================================================
    ArrayReference< T > operator = (const Array< T >& v)
    {
        Array< T >::operator = (v);
        return *this;
    }


    //========================================================================
    ArrayReference< T > copyReference(Array< T >& v)
    {
        this->d    = v.d;
        this->e    = v.e;
        this->nd   = v.nd;
        this->ne   = v.ne;
        this->stat = true;
        this->eEnd = this->e + this->ne;
        return *this;
    }

    //========================================================================
    ArrayReference< T > copyReference(ArrayReference< T > v)
    {
        this->d    = v.d;
        this->e    = v.e;
        this->nd   = v.nd;
        this->ne   = v.ne;
        this->stat = true;
        this->eEnd = this->e + this->ne;
        return *this;
    }
};

//========================================================================
template< class T >
inline ArrayReference< T > Array< T >::operator [ ](unsigned i)
{
    SIZE_CHECK(nd > 0)
    RANGE_CHECK(i < d[ 0 ])
    return ArrayReference< T > (d + 1, e + i*(ne / d[ 0 ]), nd - 1, ne / d[ 0 ]);
}


//========================================================================
template< class T >
inline const ArrayReference< T > Array< T >::operator [ ](unsigned i) const
{
    SIZE_CHECK(nd > 0)
    RANGE_CHECK(i < d[ 0 ])
    return ArrayReference< T > (d + 1, e + i*(ne / d[ 0 ]), nd - 1, ne / d[ 0 ]);
}


#endif //__ARRAY_H