---

TensorArray.H

Go to the documentation of this file.

#ifndef _TENSORARRAY_H
#define _TENSORARRAY_H


#include <vector>
#include <cstdarg>
#include <iterator>

#include <TensorArrayRowIterator.H>
#include <TensorArrayIndexIterator.H>
#include <TensorArraySliceIterator.H>
#include <arglst.H>
#include <xmath.H>
#include <ListAssignment.H>
#include <tensor_array_typedef.H>
#include <Transform.H>

using namespace std;

namespace xchen
{
  //   template<typename T, int dims>
  //   class TensorArraySliceIterator;
  
  template<typename T, int dims>
  class TensorArray
  {
      def_type_name_for_tensor_array(T, dims);

      friend class TensorArray<T,dims-1>;

  public:
      TensorArray();                                   
      explicit TensorArray(int sz, ...);               
      explicit TensorArray(int* sz);                   
      TensorArray(MyType const&);                      
      TensorArray(MyType const&, int d, int s, int e); 
      TensorArray(MyType const&, int*start, int*end);  

      TensorArray(TensorArray<T,dims+1> const& rhs, int d, int s); 


      ~TensorArray()                                   { delete [] data; delete idx_itr; }

      TensorArray& Transform(xchen::Transform const&); 
      
      void Clear();                                    
      bool Empty() const                               { return find(sz, sz+dims, 0) != sz+dims; }

      void PopFront(int i=0);                         
      void PopBack(int i=0);                           

      void Insert(linear_iterator itr, T const& val); 
      void Insert(int dir, int before);                
      
      void SetN(int N);                                

      void InsertMiddlePoints(bool average=true);      

      void Averaging();                                

      TensorArray& operator=(const TensorArray& rhs);  
      
      int GetSize(int dir=0) const                     { return sz[dir]; }  
      void GetSize(int* Sz) const                      { copy(sz, sz+dims, Sz); } 
      int GetTotalSize() const                         { return totalSz; }  
      
      void Reserve(int newCapacity,...);               
      void Reserve(int* newCapacity);                  
      void Resize(int newSz, ...);                     
      void Resize(int* newSz);                         
      void ResizeOneDirection(int i, int new_sz);      
      void CopyExtend(int left, ...);                  
      void CopyExtend(int* rightleft);                 
      void CopyExtendAtDirection(int i, 
                                 int left, int right); 
      void ExtendAtDirection(int i, 
                             int left, int right);     
      void Extend(int left, ...);                      
      void Extend(int* rightleft);                     

      slc_iterator SliceIterate(int dir=0, int start=0, int end=0)        { return slc_iterator(this, dir, start, end); }
      
      linear_iterator LinearIterate(int dir, int start, int len);         
      linear_iterator LinearIterate(int*start=0, int*end=0)               { return linear_iterator(this, start, end); }  
      linear_iterator Begin()                                             { return LinearIterate().Begin(); }            
      linear_iterator End()                                               { return LinearIterate().End(); }              

      linear_riterator LinearReverseIterate(int dir, int start, int len); 
      linear_riterator LinearReverseIterate(int*start=0, int*end=0)       { return linear_riterator(this, start, end); } 
      linear_riterator Rbegin()                                           { return LinearReverseIterate().Begin(); }     
      linear_riterator Rend()                                             { return LinearReverseIterate().End(); }       

      row_iterator RowIterate(int row_d,int subarray_d,int start,int len);
      row_iterator RowIterate(int d=0, int* start = 0, int* end=0)        { return row_iterator(d, this, start, end); }  
      row_riterator RowReverseIterate(int r_d,int sub_d,int start,int L); 
      row_riterator RowReverseIterate(int d=0, int* start=0, int* end=0)  { return row_riterator(d, this, start, end); } 

      idx_op_ret_t operator[](int i);                                     
      T& operator[](Vector<int,dims>const &i);                            
      ListAssignmentCheckSzRt<T,linear_iterator> operator=(T);            

      T const& GetElement(int* idx) const;                                
      T const& GetElement(int i, ...) const;                              
      
      friend ostream& operator<< <>(ostream& os, const TensorArray&);
      void PrintShape() const;                                            

  private:
      T* data;
      int sz[dims], totalSz;                                              // acutal elements in each dimension and total elements.
      int capacity[dims];                                                 // capacity in each dimension.
      int n;                                                              // reserved space between eles(rows) for subdivision.
      int stride[dims+1];                                                 // stride[i]: ptr_diff between ele[same]..[i]..[Same] and ele[same]..[i+1]..[Same]
      idx_itr_t* idx_itr;

  private:
      void compute_stride_totalSz();                                      // compute stride.
      void allocate()                                                     { delete [] data; data = stride[dims]? new T[stride[dims]] : 0; }
      void new_idx_itr()                                                  { delete idx_itr; idx_itr = new idx_itr_t(stride, sz); }

      void reset_after_resize();                                          // change capcaity and reallocate if necessary.
      void init_after_resize();                                           // init after only sz set: set cap = sz and allocate.
      void init_after_reshape();                                          // init after both sz and cap set.

      void initElements(const TensorArray& rhs);                          // Initialize the maximal common subarray as rhs.
      
      bool generate_even_idx_in_i_plus_dimension(int i=-1, int* start=0); // For InsertMiddlePoints(). Default call () to generate idx into start.
      bool generate_idx_of_boundary_element(int* idx=0, int*s=0,int*e=0); // For extend(). Pass paras to init; default call () to generate idx.
      void copy_extend(int* idx, int* start, int* end);                   // Copy-extend ele[idx].
      void extend(bool cpy, int* leftright);  
  };



  // Just for correct compilation.
  template<typename T>
  class TensorArray<T,0>
  {
  public:
      TensorArray(TensorArray<T,1> const&, int, int) { }
      friend  ostream&  operator<< (ostream& os, const TensorArray&) { return os; }
  };

  struct IdxRange 
  {
      IdxRange(int _d, int _s, int _e) : d(_d), s(_s), e(_e) { }
      int d, s, e;  // index range: [s, e) along direction d.
  };



  template<typename T, int dims>
  inline TensorArray<T,dims> :: TensorArray() : data(0), n(0), idx_itr(0)
  {
    bzero(sz, sizeof(int) * dims);
    bzero(capacity, sizeof(int) * dims);
    compute_stride_totalSz();
  }
  template<typename T, int dims>
  inline TensorArray<T,dims> :: TensorArray(int arg1, ...) : data(0), n(0), idx_itr(0)
  {
    get_reverse_ints_va(sz, dims);
    init_after_resize();
  }

  template<typename T, int dims>
  inline TensorArray<T,dims> :: TensorArray(int* Sz) : data(0), n(0), idx_itr(0)
  {
    copy(Sz, Sz+dims, sz);
    init_after_resize();
  }
  template<typename T, int dims>
  inline TensorArray<T,dims> :: TensorArray(MyType const& rhs) : data(0), n(0), idx_itr(0)
  {
    copy(&rhs.sz[0], &rhs.sz[0]+dims, sz);
    copy(&rhs.capacity[0], &rhs.capacity[0]+dims, capacity);
    init_after_reshape();

    TensorArray& _rhs = (TensorArray&) rhs;
    std::copy(_rhs.Begin(), _rhs.End(), Begin());
  }

  template<typename T, int dims>
  inline TensorArray<T,dims> :: TensorArray(MyType const& rhs, int d, int s, int e) : data(0), n(0), idx_itr(0)
  {
    copy(&rhs.sz[0], &rhs.sz[0]+dims, sz);
    sz[d] = e - s;

    init_after_resize();

    TensorArray& _rhs = (TensorArray&)rhs;
    linear_iterator itr = _rhs.LinearIterate(d, s, e);
    std::copy(itr.Begin(), itr.End(), Begin());
  }
  
  
  template<typename T, int dims>
  inline TensorArray<T,dims> :: TensorArray(TensorArray<T,dims+1> const& rhs, int d, int s) : data(0), n(0), idx_itr(0)
  {
    copy(&rhs.sz[0], &rhs.sz[0]+d, sz);
    copy(&rhs.sz[0]+d+1, &rhs.sz[0]+dims+1, sz);

    init_after_resize();

    SuperType& _rhs = (SuperType&)rhs;
    typename SuperType::linear_iterator itr = _rhs.LinearIterate(d, s, 1);
    std::copy(itr.Begin(), itr.End(), Begin());
  }
  
  
  template<typename T, int dims>
  inline TensorArray<T,dims> :: TensorArray(MyType const& rhs, int* s, int* e) : data(0), n(0), idx_itr(0)
  {
    transform(s, s+dims, e, sz, minus<int>());
      
    init_after_resize();

    TensorArray& _rhs = (TensorArray&)rhs;
    linear_iterator itr = _rhs.LinearIterate(s, e);
    std::copy(itr.Begin(), itr.End(), Begin());
  }
  
  
  template<typename T, int dims>
  inline TensorArray<T,dims>& TensorArray<T,dims>::Transform(xchen::Transform const& trans)
  {
    for(linear_iterator itr = Begin(); itr!=End(); ++itr)
    {
      *itr = trans(*itr);
    }
    return *this;
  }
  


  template<typename T, int dims>
  inline void TensorArray<T,dims> :: Clear()
  {
    delete [] data; delete idx_itr;

    data = 0; idx_itr = 0; n = 0;
    
    bzero(sz, sizeof(int) * dims);
    bzero(capacity, sizeof(int) * dims);
    compute_stride_totalSz();
  }

  template<typename T, int dims>
  TensorArray<T,dims>& TensorArray<T,dims> :: operator=(const TensorArray& rhs)
  {
    if(this != &rhs)
    {
      copy(&rhs.sz[0], &rhs.sz[0]+dims, sz);

      reset_after_resize();
      
      TensorArray& _rhs = (TensorArray&)rhs;
      std::copy(_rhs.Begin(), _rhs.End(), Begin());
    }
    return *this;
  }


  template<typename T, int dims>
  ListAssignmentCheckSzRt<T, typename TensorArray<T,dims> :: linear_iterator> TensorArray<T,dims> ::  operator=(T x)
  {
    data[0] = x;
    return ListAssignmentCheckSzRt<T, linear_iterator>(Begin()+1, totalSz-1);
  }

  template<typename T, int dims>
  void TensorArray<T,dims> :: initElements(const TensorArray& rhs)
  {
    if(this != &rhs)
    {
      int end[dims];
      for(int i=0; i<dims; i++)
        end[i] = std::min(sz[i], rhs.sz[i]);
        
      linear_iterator itr1 = ((TensorArray&)rhs).LinearIterate(0, end);
      linear_iterator itr2 = LinearIterate(0, end);
      std::copy(itr1.Begin(), itr1.End(), itr2.Begin());
    }
  }


  template<typename T, int dims>
  inline void TensorArray<T,dims> :: Reserve(int arg1, ...)
  {
    int new_cap[dims];
    get_reverse_ints_va(new_cap, dims);
    Reserve(new_cap);
  }
  template<typename T, int dims>
  void TensorArray<T,dims> :: Reserve(int* _new_cap)
  {
    bool need_alloc = false;
    int new_cap[dims];
    copy(_new_cap, _new_cap+dims, new_cap);
    
    for(int i=0; i<dims; i++)
    {
      if(capacity[i] < new_cap[i])
        need_alloc = true;
      else
        new_cap[i] = capacity[i];
    }

    if(need_alloc)
    {
      TensorArray<T,dims> old_ar(*this);

      copy(new_cap, new_cap+dims, capacity);
      init_after_reshape();

      copy(old_ar.Begin(), old_ar.End(), Begin());
    }
  }

  template<typename T, int dims>
  inline void TensorArray<T,dims> :: reset_after_resize()
  {
    for(int i=0; i<dims; i++)
      if(capacity[i] < sz[i])
        capacity[i] = sz[i];

    int old_total_memsz = stride[dims];
    compute_stride_totalSz();
    new_idx_itr();

    if( !data || old_total_memsz < stride[dims] )
      allocate();
  }
  template<typename T, int dims>
  inline void TensorArray<T,dims> :: init_after_resize()
  {
    copy(sz, sz+dims, capacity);
    compute_stride_totalSz();
    new_idx_itr();
    allocate();
  }
  
  template<typename T, int dims>
  inline void TensorArray<T,dims> :: init_after_reshape()
  {
    compute_stride_totalSz();
    new_idx_itr();
    allocate();
  }
  

  template<typename T, int dims>
  inline void TensorArray<T,dims> :: Resize(int arg1, ...)
  {
    int new_sz[dims];
    get_reverse_ints_va(new_sz, dims);
    Resize(new_sz);
  }
  template<typename T, int dims>
  void TensorArray<T,dims> :: Resize(int* new_sz)
  {
    TensorArray<T,dims>* old_ar = data? new TensorArray<T,dims>(*this) : 0;

    copy(new_sz, new_sz+dims, sz);
    init_after_resize();
    
    if(old_ar)
    {
      initElements(*old_ar);
      delete old_ar;
    }
    
  }

  template<typename T, int dims>
  void TensorArray<T,dims> :: ResizeOneDirection(int i, int new_sz)
  {
    if(new_sz == sz[i])
      return;
    
    TensorArray<T,dims>* old_ar = data? new TensorArray<T,dims>(*this) : 0;

    bool need_alloc = false;
    int less_sz[dims];
    GetSize(less_sz);

    if(new_sz < sz[i]) 
      less_sz[i] = new_sz;
    else 
    {
      if(new_sz > capacity[i])
      {
        capacity[i] = new_sz;
        need_alloc = true;
      }
    }
    sz[i] = new_sz;
    
    compute_stride_totalSz();
    new_idx_itr();

    if(need_alloc)
      allocate();

    if(old_ar)
    {
      linear_iterator new_itr = LinearIterate(0, less_sz);
      linear_iterator old_itr = old_ar->LinearIterate(0, less_sz);
      std::copy(old_itr.Begin(), old_itr.End(), new_itr.Begin());
      delete old_ar;
    }
    
  }

  template<typename T, int dims>
  void TensorArray<T,dims> :: Insert(linear_iterator itr, T const& val) 
  {
    T* cpy = new T[totalSz];
    T *bgn = cpy, *inst, *end = cpy+totalSz;

    linear_iterator it = Begin();
    for(; it != itr; ++it, ++cpy) *cpy = *it;

    inst = cpy;

    for(; it != End(); ++it, ++cpy) *cpy = *it;

    ResizeOneDirection(dims-1, sz[dims-1]+1);
    
    cpy = bgn;

    for(it = Begin(); bgn != inst; ++bgn, ++it) *it = *bgn;
    
    *it = val; ++it;
    
    for(; bgn != end; ++bgn, ++it) *it = *bgn;

    delete [] cpy;
  }
  

  template<typename T, int dims>
  void TensorArray<T,dims> :: Insert(int dir, int before)
  {
    int len = sz[dir];
    TensorArray<T,dims> cpy(*this);
    ResizeOneDirection(dir, len+1);
    array_copy(dir,before,  cpy,0,  *this,0);
    array_copy(dir,len-before,  cpy,before,  *this,before+1);
  }
  
  
  template<typename T, int dims>
  inline typename TensorArray<T,dims> :: linear_iterator TensorArray<T,dims> :: LinearIterate(int dir, int start_idx, int len)
  {
    int start[dims], end[dims];
    bzero(start, sizeof(int) * dims); GetSize(end);
    start[dir] = start_idx, end[dir] = start_idx + len;
    return LinearIterate(start, end);
  }
  template<typename T, int dims>
  inline typename TensorArray<T,dims> :: linear_riterator TensorArray<T,dims> :: LinearReverseIterate(int dir, int start_idx, int len)
  {
    int start[dims], end[dims];
    bzero(start, sizeof(int) * dims); GetSize(end);
    start[dir] = start_idx, end[dir] = start_idx + len;
    return LinearReverseIterate(start, end);
  }
  template<typename T, int dims>
  inline typename TensorArray<T,dims> :: row_iterator TensorArray<T,dims> :: RowIterate(int row_dir, int subarray_dir, int start_idx, int len)
  {
    int start[dims], end[dims];
    bzero(start, sizeof(int) * dims); GetSize(end);
    start[subarray_dir] = start_idx, end[subarray_dir] = start_idx + len;
    return RowIterate(row_dir, start, end);
  }
  template<typename T, int dims>
  inline typename TensorArray<T,dims> :: row_riterator TensorArray<T,dims> :: RowReverseIterate(int row_dir, int subarray_dir, int start_idx, int len)
  {
    int start[dims], end[dims];
    bzero(start, sizeof(int) * dims); GetSize(end);
    start[subarray_dir] = start_idx, end[subarray_dir] = start_idx + len;
    return RowReverseIterate(row_dir, start, end);
  }
  

  template<typename T, int dims>
  void TensorArray<T,dims> :: PopFront(int i)
  {
    TensorArray<T,dims> cpy(*this);
    sz[i]--;
    compute_stride_totalSz();
    new_idx_itr();

    array_copy(i, sz[i], cpy, 1, *this, 0);
  }
  template<typename T, int dims>
  void TensorArray<T,dims> :: PopBack(int i)
  {
    sz[i]--;
    compute_stride_totalSz();
    new_idx_itr();
  }

  template<typename T, int dims>
  inline void TensorArray<T,dims> :: Extend(int arg1, ...)
  {
    int leftright[2*dims];
    get_reverse_ints_va(leftright, 2*dims);
    Extend(leftright);
  }
  template<typename T, int dims>
  inline void TensorArray<T,dims> :: Extend(int* rightleft)
  {
    extend(false, rightleft);
  }
  template<typename T, int dims>
  inline void TensorArray<T,dims> :: CopyExtend(int arg1, ...)
  {
    int rightleft[2*dims];
    get_reverse_ints_va(rightleft, 2*dims);
    CopyExtend(rightleft);
  }
  template<typename T, int dims>
  inline void TensorArray<T,dims> :: CopyExtend(int* rightleft)
  {
    extend(true, rightleft);
  }
  template<typename T, int dims>
  inline void TensorArray<T,dims> :: CopyExtendAtDirection(int i, int left, int right)
  {
    int rightleft[2*dims];
    bzero(rightleft, sizeof(int) * dims * 2);
    rightleft[2*i] = right, rightleft[2*i+1] = left;
    CopyExtend(rightleft);
  }
  template<typename T, int dims>
  inline void TensorArray<T,dims> :: ExtendAtDirection(int i, int left, int right)
  {
    int rightleft[2*dims];
    bzero(rightleft, sizeof(int) * dims * 2);
    rightleft[2*i] = right, rightleft[2*i+1] = left;
    Extend(rightleft);
  }
    

  template<typename T, int dims>
  void TensorArray<T,dims> :: extend(bool cpy_extend, int* rightleft)
  {
    TensorArray<T,dims> cpy(*this);
    int newsz[dims];
    int start[dims], end[dims];                       // old array as the subarray in the extended array.
    copy(sz, sz+dims, newsz);
    for(int i=0; i<dims; i++)
    {
      newsz[i] += rightleft[2*i] + rightleft[2*i+1];
      start[i] = rightleft[2*i+1];
      end[i] = start[i] + sz[i];
    }
    Resize(newsz);
    linear_riterator itr = LinearReverseIterate(start, end);
    copy(cpy.Rbegin(), cpy.Rend(), itr.Begin());

    if(cpy_extend) 
    {
      int idx[dims];
      generate_idx_of_boundary_element(idx, start, end);
      do 
      {
        copy_extend(idx, start, end);
      } while( generate_idx_of_boundary_element() );
    }
  }

  /* iterate all boundary indices */
  template<typename T, int dims>
  bool TensorArray<T,dims> :: generate_idx_of_boundary_element(int* idx, int* start, int* end)
  {
    const int at_interior = 0, at_start = 1, at_end = 2;  // where is the index[] value? at_end actually means right before end. 
    static int *index, s[dims], e[dims], at[dims], len[dims];
    
    if(idx && start && end)
    {
      index = idx;
      copy(start, start+dims, s);
      copy(end, end+dims, e);
      transform(e, e+dims, s, len, minus<int>());

      /* start with all 0s, i.e. all at_interiro except mandatory at_start or at_start at rightest position. */      
      for(int i=0; i<dims; i++)
      {
        if(len[i]>2)
        {
          at[i] = at_interior; index[i] = s[i] + 1;
        }
        else
        {
          at[i] = at_start; index[i] = s[i];      // no at_interior for this position.
        }
      }
      if(find(at, at+dims, at_start) == at+dims)
      {
        at[0] = at_start;  index[0] = s[0];
      }
      
      return true;
    }

    while(1) 
    {
      int i; 
      for(i=0; i<dims; i++)    // iterate all interior points with regarding to the current set of fixed end indices.
      {
        if(at[i] == at_interior)
        {
          if(++index[i] < e[i]-1)
            return true;
          index[i] = s[i] + 1;
        }
      }
      if(i==dims)             // generate next set of fixed end indices.
      {
        for(i=0; i<dims; i++)
        {
          if(len[i]==1) continue;
          ++at[i];
          if(at[i] == at_start)       { index[i] = s[i]; break; }
          else if(at[i] == at_end)    { index[i] = e[i]-1; break; }
          else  // carry occurs. 
            if(len[i]>2)              { at[i] = at_interior; index[i] = s[i] + 1; } 
            else                      { at[i] = at_start; index[i] = s[i]; }
        }
      }
      return i!=dims;
    }
  }


  
  template<typename T, int dims>
  void TensorArray<T,dims> :: copy_extend(int* idx, int* start, int* end)
  {
    int s[dims], e[dims];
    copy(idx, idx+dims, s);
    transform(idx, idx+dims, e, bind2nd(plus<int>(), 1)); // initialize to-be-extended-subarray as the element itself.

    bool need_ext = false;
    for(int i=0; i<dims; i++)
    {
      if(start[i] && idx[i] == start[i])
      {
        s[i] = 0; need_ext = true;
      }
      if(idx[i] == end[i]-1 && end[i] != sz[i])
      {
        e[i] = sz[i]; need_ext = true;
      }
    }

    if(need_ext)
    {
      linear_iterator itr = LinearIterate(s,e);
      fill(itr.Begin(), itr.End(), GetElement(idx));
    }
  }
  
      
    
  // improvment ....
  template<typename T, int dims>
  void TensorArray<T,dims> :: InsertMiddlePoints(bool average)
  {
    if(!n) SetN(1);

    assert( find(sz, sz+dims, 1) == sz+dims );
    
    for(int i=0; i<dims; i++)
    {
      sz[i] = 2*sz[i] - 1;
      capacity[i] = 2 * capacity[i];
    }
    n--;
    compute_stride_totalSz();
    new_idx_itr();

    if(!average) return;
    
    int subStart[dims], subEnd[dims];
    for(int i=0; i<dims; i++)
    {
      subStart[i] = 0;
      subEnd[i] = sz[i];
      generate_even_idx_in_i_plus_dimension(i, subStart);   // initialization.
      while(generate_even_idx_in_i_plus_dimension())
      {
        for(int j=i+1; j<dims; j++)
          subEnd[j] = subStart[j] + 1;
        
        for( row_iterator r_itr = RowIterate(i, subStart, subEnd); !r_itr.PastEnd(); ++r_itr)
        {
          element_iterator e_itr_left = r_itr.Begin();
          element_iterator e_itr_right = r_itr.Begin() + 2;
          element_iterator e_itr_middle = r_itr.Begin() + 1;
          for(; e_itr_right < r_itr.End(); e_itr_left += 2, e_itr_right += 2, e_itr_middle += 2)
          {
            *e_itr_middle = (*e_itr_left + *e_itr_right) / 2;
          }
        }
      }
    }
  }

  template<typename T, int dims>
  bool TensorArray<T,dims> :: generate_even_idx_in_i_plus_dimension(int ii, int * subStart)
  {
    static int i, *sub_start;
    static bool first = false;
    
    if(subStart)
    {
      i = ii+1;
      sub_start = subStart;
      for(int j=i; j<dims; j++)
        sub_start[j] = j==i? -2 : 0;
      first = true;
      return true;
    }
    else if(i >= dims)
      return first? (first=false, true) : false;
    else
    {
      for(int j=i; j<dims; j++)
      {
        if( (sub_start[j] += 2) < sz[j] )
          return true;
        else
          sub_start[j] = 0;
      }

      return false;
    }
  }
  
        
  template<typename T, int dims>
  inline void TensorArray<T,dims> :: SetN(int N)  
  {
    if(n!=N) 
    {
      TensorArray<T,dims> cpy(*this);
      n = N;
      compute_stride_totalSz();
      allocate();
      new_idx_itr();
    
      std::copy(cpy.Begin(), cpy.End(), Begin());
    }
  }

  /* Average each hyper cube into its 'lowerleft' corner.*/
  template<typename T, int dims>
  void TensorArray<T,dims> :: Averaging()
  {
    int eles_per_cell = int_pow(2, dims);
    linear_iterator itr[eles_per_cell];
    int start[dims], len[dims], end[dims];           // total eles_per_cell subarrays.

    bzero(start, sizeof(int) * dims);
    GetSize(len);
    transform(len, len+dims, len, bind2nd(minus<int>(), 1));

    
    /*
     * Generate indices of each vertex of the first hyper-cube via binaray number representation. 
     * The linear iterator for the vertex is generated based on that.
     */
    for(int i=0; i<eles_per_cell; i++)  
    {
      for(int j = 0; i && !(start[j] = 1 - start[j]) && j<dims; j++);
      transform(start, start+dims, len, end, plus<int>());
      itr[i] = LinearIterate(start, end);
    }
    
    /* Average while iterating the corresponding vertices of the sequential hypercubes. */
    while( !itr[0].PastEnd() )
    {
      for(int i=1; i<eles_per_cell; i++)
      {
        *itr[0] += *itr[i];          
      }
      *itr[0] /= eles_per_cell;
      for(int i=0; i<eles_per_cell; i++)
        ++itr[i];
    }
    GetSize(len);
    transform(len, len+dims, len, bind2nd(minus<int>(), 1));
    Resize(len);
  }
  

  

  template<typename T, int dims>
  inline void TensorArray<T,dims> :: compute_stride_totalSz()
  {
    int strd = (1 << n);
    stride[0] = strd;
    for(int i=1; i<=dims; i++)
      stride[i] = stride[i-1] * capacity[i-1]  * strd;
    totalSz = sz[0];
    for(int i=1; i<dims; i++)
      totalSz *= sz[i];
  }


  template<typename T, int dims>
  inline T& TensorArray<T,dims> :: operator[](Vector<int,dims>const& idx)
  {
    T* ptr = data;
    for(int i=dims-1; i>=0; i--)
    {
      ptr += stride[i] * idx[i];
    }
    return *ptr;
  }
  

  template<typename T, int dims>
  inline typename TensorArrayIndexOperatorTraits<T,dims>::ReturnType TensorArray<T,dims> :: operator[](int i)          
  { 
    assert(i<sz[dims-1]);
    if(dims>1) 
    {
      idx_itr->setPtr(data); 
      return (*idx_itr)[i];
    }
    return  idx_op_ret_t(data[i]); 
  }

  template<typename T, int dims>
  T const& TensorArray<T,dims> :: GetElement(int arg1, ...) const
  {
    int idx[dims];
    get_reverse_ints_va(idx, dims);
    return GetElement(idx);
  }
  template<typename T, int dims>
  T const& TensorArray<T,dims> :: GetElement(int* idx) const
  {
    T* ptr = data;
    for(int i=0; i<dims; i++)
      ptr += idx[i] * stride[i];
    return *ptr;
  }
  

  template<typename T, int dims>
  void TensorArray<T,dims> :: PrintShape() const
  {
    cout << "\nTensorArray of \nSize:\t\t";
    for(int i=dims-1; i>=0; i--)  cout << "["<< sz[i] << "]";
    cout << "\nCapacity:\t";
    for(int i=dims-1; i>=0; i--)  cout << "["<< capacity[i] << "]";
    cout << "\nReserve space:\t2^" << n << "-1\n";
  }


  template<typename T, int dims>
  ostream&  operator<< (ostream& os, const TensorArray<T, dims>& _ar)
  {
    TensorArray<T, dims>& ar = (TensorArray<T, dims>&) _ar;
    if(dims==1)
      std::copy(ar.Begin(), ar.End(), ostream_iterator<T>(os, " "));
    else 
    {
      for(int start = 0; start < ar.GetSize(dims-1); ++start)
      {
        TensorArray<T, dims-1> sub_array(ar, dims-1, start);
        os << sub_array << endl;
      }
    }
    return os ;
  }



  template<typename T, int dims1, int dims2>
  inline void array_copy(int sub_array_at_this_dir, int len, 
                         TensorArray<T,dims1>const& from, int from_start, 
                         TensorArray<T,dims2>& to, int to_start)
  {
    int i = sub_array_at_this_dir;
    int fromStart[dims1], fromEnd[dims1], toStart[dims2], toEnd[dims2];

    bzero(fromStart, sizeof(int) * dims1);
    from.GetSize(fromEnd); 
    fromStart[i] = from_start, fromEnd[i] = from_start + len;

    bzero(toStart, sizeof(int) * dims2);
    to.GetSize(toEnd); 
    toStart[i] = to_start; toEnd[i] = to_start + len;

    typename TensorArray<T,dims1> :: linear_iterator  
      from_itr = ( (TensorArray<T,dims1>&)(from) ).LinearIterate(fromStart, fromEnd);
    
    typename TensorArray<T,dims2> :: linear_iterator  
      to_itr = to.LinearIterate(toStart, toEnd);

    copy(from_itr.Begin(), from_itr.End(), to_itr.Begin());
  }
  


  template<typename T, int dims>
  inline TensorArray<T,dims> operator*(Transform const& trans, TensorArray<T,dims>const& ar)
  {
    TensorArray<T,dims> ret(ar);
    ret.Transform(trans);
    return ret;
  }
  


}//end namespace xchen


#endif

---