//
//  Dim.h
//
//  Qdos2 - Dim Class used for defining dimensionality of a Matrix Object
//
//  Janet Barnes 
//  Sept 1999
//
#ifndef QDOS_DIM_H
#define QDOS_DIM_H

#include <vector>
#include <valarray>
#include "Exception.h"
#include "Viostream.h"
#include "Var.h"
#include <numeric>
using namespace std;

namespace QDOS {

typedef enum DimStorageOption_e {ROW_MAJOR_ORDER,COL_MAJOR_ORDER} \
                DimStorageOption;

// ++++++++++++++++++++++++++++++++++++++++++++++++++ misc. functions


    inline valarray<size_t> computestrides(size_t *dm,size_t rk,\
                                    DimStorageOption order=ROW_MAJOR_ORDER)
    // calculate strides betweeen dimensions..
    {
        valarray<size_t> st(1,rk);
        if (rk>0)
        {
            if (order==ROW_MAJOR_ORDER) // row-major order (default)
                for (int i=rk-1;i>0;i--){st[i-1]=st[i]*dm[i];}
            else    //col-major order
                for (int i=0;i<rk-1;i++){st[i+1]=st[i]*dm[i];}
        }
        return st;
    } 
    
    inline valarray<size_t> computestrides(const valarray<size_t>& dm,\
                                    DimStorageOption order=ROW_MAJOR_ORDER)
    // calculate strides betweeen dimensions..
    {
        size_t rk=dm.size();
        valarray<size_t> st(1,rk);
        if (rk>0)
        {
            if (order==ROW_MAJOR_ORDER) // row-major order (default)
                for (int i=rk-1;i>0;i--){st[i-1]=st[i]*dm[i];}
            else    //col-major order
                for (int i=0;i<rk-1;i++){st[i+1]=st[i]*dm[i];}    
        }
        return st;
    }
    
    inline valarray<size_t> setdims(size_t *dms,size_t ndims)
    {
    // set dimension info
       return valarray<size_t>(dms,ndims);
    } 

    inline size_t product(valarray<size_t> v)
    { return std::accumulate(&v[0],&v[v.size()],size_t(1),multiplies<size_t>());}
	    
// ----------- End of misc functions ------------------- 

class Dim;
Vostream& operator<<(Vostream&,const Dim&);
Vistream& operator>>(Vistream&, Dim&);

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ class Dim

// Dimension info for a Matrix object
//-------------------------------------------------------------------------    
class Dim:public VarTarget
{  
  friend Vostream& operator<<(Vostream&,const Dim&);
  friend Vistream& operator>>(Vistream&, Dim&);
private:
    gslice _gsl;                    // start index, dims and strides
    size_t _size;                    // total no. of elts
    size_t _rank;                    // rank
    DimStorageOption _order;         // whether in row major or col major order
          
public:
    // --- Constructor functions
    Dim(DimStorageOption order=ROW_MAJOR_ORDER):
        _rank(0),_order(order),_gsl(initgslice()),_size(1){;} // rank 0
  
     
    Dim(size_t d1,DimStorageOption order=ROW_MAJOR_ORDER): 
        _rank(1),_order(order),_gsl(initgslice(d1,order)),
        _size(product(_gsl.size()))
        {;}// rank 1
   
    Dim(size_t d1,size_t d2,DimStorageOption order=ROW_MAJOR_ORDER):
        _rank(2),_order(order),_gsl(initgslice(d1,d2,order)),
        _size(product(_gsl.size()))
       {;} // rank 2 
       
    Dim(size_t d1,size_t d2,size_t d3,DimStorageOption order=ROW_MAJOR_ORDER):
        _rank(3),_order(order),_gsl(initgslice(d1,d2,d3,order)),
        _size(product(_gsl.size()))
        {;} // rank 3
 
    Dim(size_t *dms,size_t nrank,DimStorageOption order=ROW_MAJOR_ORDER):
        _rank(nrank),_order(order),_gsl(initgslice(dms,nrank,order))
    {
        _size=product(_gsl.size());
    }
        
    Dim(const valarray<size_t>& dms, DimStorageOption order=ROW_MAJOR_ORDER);

    Dim(const vector<int>& dms, DimStorageOption order=ROW_MAJOR_ORDER);
    
    Dim(size_t start, const valarray<size_t>& dms, \
        const valarray<size_t>& strides, DimStorageOption order=ROW_MAJOR_ORDER):\
        _order(order),_rank(dms.size())
    {
        if(_rank==0)
        {
            valarray<size_t> v((size_t)1,1);
            _gsl=gslice(start,v,v);
        }
        else 
            _gsl=gslice(start,dms,strides);
        _size=product(_gsl.size());
    }    
       
    Dim(const gslice& g,DimStorageOption order=ROW_MAJOR_ORDER):
        _order(order)
    {
       valarray<size_t> sz(g.size());
       _rank=0;
       valarray<bool> msk(true,sz.size());
       for(int i=0;i<sz.size();i++)
       {
          msk[i]=(sz[i]!=0);
          if(msk[i]) _rank++;
       }
    // construct from gslice info
        if (_rank==0)
        {
            valarray<size_t> v((size_t)1,1);
            _gsl=gslice(g.start(),v,v);
        }
        else
        {
           _gsl=gslice(g.start(),valarray<size_t>(g.size()[msk]),\
             valarray<size_t>(g.stride()[msk]));
        }
       _size=product(_gsl.size());
    }    
           
    Dim(const Dim& d):_size(d._size),_rank(d._rank),_order(d._order),
        _gsl(d._gsl)    //_gsl(initgslice(d))
       {;} // copy constructor


    Dim& operator=(const Dim& d)
    {  _size=d._size;
       _rank=d._rank;
       _order=d._order;
       _gsl=d._gsl; 
     } // assignment operator

    // --- Operators for indexing []
    size_t operator[](const size_t& i) const
    { return _gsl.size()[i];}
            
    // --- Logical operators
    bool operator==(const Dim& d) const
    {
         if (rank()!=d.rank())
            return false;
         valarray<bool> msk(dims()!=d.dims());
         valarray<size_t> v(dims()[msk]);
         if (v.size()!=0)
             return false;
         else
            return true;
    }

    ~Dim() 
    {;}
    
    bool operator !=(const Dim& d) const
    {
        return !((*this)==d);
    }
        
    // --- Miscellaneous functions
    
    static inline gslice initgslice()
    {
        size_t one=1;
        valarray<size_t> v(one,1);
        return gslice(0,v,v);
    }
    
    static inline gslice initgslice(size_t d1,DimStorageOption order=ROW_MAJOR_ORDER)
    {
         size_t dms[1]={d1};
         valarray<size_t> vv(setdims(dms,1));
         return gslice(0,setdims(dms,1),computestrides(dms,1,order));
    }
    
    static inline gslice initgslice(size_t d1,size_t d2,
                      DimStorageOption order=ROW_MAJOR_ORDER)
    {
         size_t dms[2]={d1,d2};
         return gslice(0,setdims(dms,2),computestrides(dms,2,order));
    }
    
    static inline gslice initgslice(size_t d1,size_t d2,size_t d3,
                      DimStorageOption order=ROW_MAJOR_ORDER)
    {
         size_t dms[3]={d1,d2,d3};
         return gslice(0,setdims(dms,3),computestrides(dms,3,order));
    }
    
    static inline gslice initgslice(size_t* dms,size_t ndims,
                      DimStorageOption order=ROW_MAJOR_ORDER)
    {
         if(ndims==0)
             return initgslice();
         return gslice(0,setdims(dms,ndims),computestrides(dms,ndims,order));
    }
  
    static inline gslice initgslice(const valarray<size_t>& dms,
                      DimStorageOption order=ROW_MAJOR_ORDER)
    {
        if(dms.size()==0)
            return initgslice();
        return gslice(0,dms,computestrides(dms,order));
    }
    
    bool is_contiguous()
    {
    // are indices as specified by gslice _gsl contiguous?
        bool res=1;
        
        valarray<size_t> sz(dims());
        valarray<size_t> st(computestrides(sz,_order));
        
        for (size_t i=0;i<rank();i++)
             res=res && strides()[i]==st[i];
        return res;
    }
    
    bool is_square() const
    {
         size_t i;
         bool res=1;
         
         // check that all elements of dims() are equal
         for (i=1;i<rank();i++)
             res=res && dims()[0]==dims()[i];
         return res;
    }    
    
    DimStorageOption storageorder() const{return _order;}
    valarray<size_t> dims() const
    { if (_rank==0) return valarray<size_t>();
      return _gsl.size();} // dimension info
    valarray<size_t> strides() const{return _gsl.stride();} // strides
    void get_strides(valarray<size_t>& str) const{str=_gsl.stride();}
    void get_strides(vector<int>&) const;
    void get_dims(valarray<size_t>& dm) const{dm=_gsl.size();}
    void get_dims(vector<int>&) const;
    size_t size() const{return _size;}  // total no. elts
    size_t rank() const{return _rank;}  // no. of dimensions
    size_t start() const{return _gsl.start();} // start index 
    gslice& getgslice() {return _gsl;}     // return gslice of info

    bool samesize(const Dim& d) const; // Dim Objects identical in size & rank
    void display() const; // display function
        
}; 
// ----------- End of class Dim -------------------

typedef Var<Dim> Dim_var;
inline Vostream& operator<<(Vostream& v_os ,const Dim& dm)
{ 
  v_os << dm._gsl;
  v_os << dm._size;
  v_os << dm._rank;
  v_os << (int)dm._order;
  return v_os;
}

inline Vistream& operator>>(Vistream& v_is, Dim& dm)
{
   v_is >> dm._gsl;
   v_is >> dm._size;
   v_is >> dm._rank;
   int i;
   v_is >> i;
   dm._order=(DimStorageOption)i;
   return v_is;
}

ostream& operator<<(ostream& os,const Dim& d);

} // end namespace QDOS

#endif
