data.h

/***************************************************************************
                          data.h  -  description
                             -------------------
    begin                : Fri Apr 6 2001
    copyright            : (C) 2001 by Thies Jochimsen & Michael von Mengershausen
    email                : jochimse@cns.mpg.de  mengers@cns.mpg.de
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#ifndef DATA_H
#define DATA_H

#ifndef TJUTILS_CONFIG_H
#define TJUTILS_CONFIG_H
#include <tjutils/config.h>
#endif

#include <tjutils/tjfeedback.h>

#ifdef ODIN_DEBUG
#define BZ_DEBUG
#endif




#include <blitz/blitz.h> // include 1st to overwrite mutex defines

// undef mutex defines
#ifdef BZ_MUTEX_DECLARE
#undef BZ_MUTEX_DECLARE
#endif
#ifdef BZ_MUTEX_INIT
#undef BZ_MUTEX_INIT
#endif
#ifdef BZ_MUTEX_LOCK
#undef BZ_MUTEX_LOCK
#endif
#ifdef BZ_MUTEX_UNLOCK
#undef BZ_MUTEX_UNLOCK
#endif
#ifdef BZ_MUTEX_DESTROY
#undef BZ_MUTEX_DESTROY
#endif

// do not use mutex in Blitz, even if pthread is available
// to get consistent behaviour across all platforms
#define BZ_MUTEX_DECLARE(name)
#define BZ_MUTEX_INIT(name)
#define BZ_MUTEX_LOCK(name)
#define BZ_MUTEX_UNLOCK(name)
#define BZ_MUTEX_DESTROY(name)



#include <blitz/array.h>
#include <blitz/tinyvec-et.h> // to use arithmetics of TinyVectors
using namespace blitz;


#include <tjutils/tjtypes.h>
#include <tjutils/tjprofiler.h>

#include <odinpara/protocol.h>

#include <odindata/converter.h>
#include <odindata/fileio_opts.h>


// some additional macros for complex numbers
BZ_DECLARE_FUNCTION_RET (float2real, STD_complex);
BZ_DECLARE_FUNCTION_RET (float2imag, STD_complex);
BZ_DECLARE_FUNCTION_RET (creal, float)
BZ_DECLARE_FUNCTION_RET (cimag, float)
BZ_DECLARE_FUNCTION_RET (cabs, float)
BZ_DECLARE_FUNCTION_RET (phase, float)
BZ_DECLARE_FUNCTION (logc)
BZ_DECLARE_FUNCTION (expc)
BZ_DECLARE_FUNCTION (conjc)

// Other useful macros
BZ_DECLARE_FUNCTION2(secureDivision)
BZ_DECLARE_FUNCTION (secureInv)


#ifdef STL_REPLACEMENT
BZ_DECLARE_ARRAY_ET_SCALAR_OPS(tjstd_complex)
#endif




struct FileMapHandle {
  FileMapHandle() : fd(-1), offset(0), refcount(1) {}
  int fd;
  LONGEST_INT offset; // in bytes
  int refcount;
  Mutex mutex; // protect refcount
};


template <typename T, int N_rank>
class Data : public Array<T,N_rank> {

public:


  Data(const TinyVector<int,N_rank>& dimvec, const T& val=0) : Array<T,N_rank>(dimvec), fmap(0) {(*this)=val;}


  // Variable extent constructors
  Data() : fmap(0) {}
  Data(int extent1) : Array<T,N_rank>(extent1), fmap(0)  {}
  Data(int extent1, int extent2) : Array<T,N_rank>(extent1,extent2), fmap(0)  {}
  Data(int extent1, int extent2,int extent3) : Array<T,N_rank>(extent1,extent2,extent3), fmap(0)  {}
  Data(int extent1, int extent2,int extent3,int extent4) : Array<T,N_rank>(extent1,extent2,extent3,extent4), fmap(0)  {}
  Data(int extent1, int extent2,int extent3,int extent4,int extent5) : Array<T,N_rank>(extent1,extent2,extent3,extent4,extent5), fmap(0)  {}
  Data(int extent1, int extent2,int extent3,int extent4,int extent5,int extent6) : Array<T,N_rank>(extent1,extent2,extent3,extent4,extent5,extent6), fmap(0)  {}
  Data(int extent1, int extent2,int extent3,int extent4,int extent5,int extent6,int extent7) : Array<T,N_rank>(extent1,extent2,extent3,extent4,extent5,extent6,extent7), fmap(0)  {}
  Data(int extent1, int extent2,int extent3,int extent4,int extent5,int extent6,int extent7,int extent8) : Array<T,N_rank>(extent1,extent2,extent3,extent4,extent5,extent6,extent7,extent8), fmap(0)  {}
  Data(int extent1, int extent2,int extent3,int extent4,int extent5,int extent6,int extent7,int extent8,int extent9) : Array<T,N_rank>(extent1,extent2,extent3,extent4,extent5,extent6,extent7,extent8,extent9), fmap(0)  {}


  Data(const STD_string& filename, bool readonly, const TinyVector<int,N_rank>& shape, LONGEST_INT offset=0 );

  Data(const Data<T,N_rank>& d) : fmap(0) {Data<T,N_rank>::reference(d);}

  Data(const Array<T,N_rank>& a) : Array<T,N_rank>(a), fmap(0) {}

  Data(const tjarray<tjvector<T>,T>& a) : fmap(0) { (*this)=a;}

  template<class T_expr>
  Data(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr) : Array<T,N_rank>(expr), fmap(0) {}


  ~Data();


  Data<T,N_rank>& operator = (const Data<T,N_rank>& d) {Array<T,N_rank>::operator=(d); return *this;}

  Data<T,N_rank>& operator = (const Array<T,N_rank>& a) {Array<T,N_rank>::operator=(a); return *this;}

  Data<T,N_rank>& operator = (const T& val) {Array<T,N_rank>::operator=(val); return *this;}

  template<class T_expr>
  inline Data<T,N_rank>& operator = (BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr) {
    typedef _bz_typename T_expr::T_numtype T_numtype;
    evaluate(expr, _bz_update<T_numtype, _bz_typename T_expr::T_numtype>());
    return *this;
  }

  Data<T,N_rank>& operator = (const tjarray<tjvector<T>,T>& a);


  // resolve ambiguities
  Array<T,N_rank> operator + (const Data<T,N_rank>& b) const {return Array<T,N_rank>(Array<T,N_rank>(*this)+Array<T,N_rank>(b));}
  Array<T,N_rank> operator - (const Data<T,N_rank>& b) const {return Array<T,N_rank>(Array<T,N_rank>(*this)-Array<T,N_rank>(b));}
  Array<T,N_rank> operator * (const Data<T,N_rank>& b) const {return Array<T,N_rank>(Array<T,N_rank>(*this)*Array<T,N_rank>(b));}
  Array<T,N_rank> operator / (const Data<T,N_rank>& b) const {return Array<T,N_rank>(Array<T,N_rank>(*this)/Array<T,N_rank>(b));}


  // resolve ambiguities
  Array<T,N_rank> operator + (const T& v) const {return Array<T,N_rank>(Array<T,N_rank>(*this)+v);}
  Array<T,N_rank> operator - (const T& v) const {return Array<T,N_rank>(Array<T,N_rank>(*this)-v);}
  Array<T,N_rank> operator * (const T& v) const {return Array<T,N_rank>(Array<T,N_rank>(*this)*v);}
  Array<T,N_rank> operator / (const T& v) const {return Array<T,N_rank>(Array<T,N_rank>(*this)/v);}



  int read(const STD_string& format,const STD_string& filename, LONGEST_INT offset=0);


  template<typename T2>
  int read(const STD_string& filename, LONGEST_INT offset=0);


  int write(const STD_string& format,const STD_string& filename, autoscaleOption scaleopt=autoscale) const;


  template<typename T2>
  int write(const STD_string& filename, autoscaleOption scaleopt=autoscale) const;



  int write(const STD_string& filename, fopenMode mode=overwriteMode) const;

  int write_asc_file(const STD_string& filename, const Array<T,N_rank>& pre=defaultArray, const Array<T,N_rank>& post=defaultArray) const;

  int read_asc_file(const STD_string& filename);


  int autowrite(const STD_string& filename, const FileWriteOpts& opts=FileWriteOpts(), const Protocol* prot=0) const;

  int autoread(const STD_string& filename, const FileReadOpts& opts=FileReadOpts(), Protocol* prot=0);

  bool is_filemapped() const {return fmap;}

  operator tjarray<tjvector<T>,T> () const;


  TinyVector<int, N_rank> create_index(unsigned long index) const;

  unsigned long create_linear_index(const TinyVector<int, N_rank>& indexvec) const;

  void shift(unsigned int shift_dim, int shift);

  void congrid(const TinyVector<int,N_rank>& newshape, const TinyVector<float,N_rank>* subpixel_shift=0, bool left_to_right=false);


  template<typename T2, int N_rank2> Data<T2,N_rank2>& convert_to(Data<T2,N_rank2>& dst, autoscaleOption scaleopt=autoscale) const;

  // specialization which uses reference if type/rank are same and scale/offset==0.0 to reduce memory usage
  Data<T,N_rank>& convert_to(Data<T,N_rank>& dst, autoscaleOption scaleopt=autoscale) const;


  T* c_array();


  // Reimplemented Blitz::Array functions
  void reference(const Data<T,N_rank>& d);


  // dummy objects used for default arguments
  static Array<T,N_rank> defaultArray;


 private:

  void interpolate1dim(unsigned int dim, int newsize, float subpixel_shift);

  void detach_fmap();

  FileMapHandle* fmap;

};

template <typename T, int N_rank>
Array<T,N_rank> Data<T,N_rank>::defaultArray;




#ifdef STREAM_REPLACEMENT

template<typename T, int N_rank>
inline STD_ostream& operator << (STD_ostream& s,const Array<T,N_rank>&) {return s;}

template<typename T, int N_rank>
inline STD_ostream& operator << (STD_ostream& s,const TinyVector<T,N_rank>&) { return s;}

inline STD_ostream& operator << (STD_ostream& s,const Range&) { return s;}

#endif



/*
 * Create n-dim index from linear 'index' for a given 'shape'
 */
template<int N_rank>
TinyVector<int,N_rank> index2extent(const TinyVector<int,N_rank>& shape, unsigned int index) {
  TinyVector<int,N_rank> result;
  unsigned int temp=index;
  for(int i=N_rank-1;i>=0;i--){
    result(i)=temp%shape(i);
    temp=temp/shape(i);
  }
  return result;
}




// leave this function outside of class Data, otherwise it screws up
// when using optimization
template<typename T, int N_rank, typename T2>
void convert_from_ptr(Data<T,N_rank>& dst, const T2* src, const TinyVector<int, N_rank>& shape, autoscaleOption scaleopt=autoscale) {
  Log<OdinData> odinlog("Data","convert_from_ptr");

  int dstsize=product(shape);
  int srcsize=dstsize*Converter::get_elements((T)0)/Converter::get_elements((T2)0);
  ODINLOG(odinlog,normalDebug) << "dstsize/srcsize=" << dstsize << "/" << srcsize << STD_endl;

  dst.resize(shape);

  Converter::convert_array(src, dst.c_array(), srcsize, dstsize, scaleopt);
}

// specialization in case the type is the same
template<typename T, int N_rank>
void convert_from_ptr(Data<T,N_rank>& dst, const T* src, const TinyVector<int, N_rank>& shape) {
  dst.reference(Array<T,N_rank>((T*)src, shape, duplicateData));
}



template <typename T, int N_rank>
Data<T,N_rank>::Data(const STD_string& filename, bool readonly, const TinyVector<int,N_rank>& shape, LONGEST_INT offset ) : fmap(new FileMapHandle) {
  T* ptr=(T*)filemap(filename, (LONGEST_INT)product(shape)*sizeof(T), offset, readonly, fmap->fd);
  if(ptr && (fmap->fd)>=0) { // only negative file descriptor indicates filemap failure
    Array<T,N_rank>::reference(Array<T,N_rank>(ptr, shape, neverDeleteData));
    fmap->offset=offset;
  } else {
    delete fmap;
    fmap=0;
  }
}


template <typename T, int N_rank>
void Data<T,N_rank>::detach_fmap() {
  Log<OdinData> odinlog("Data","detach_fmap");
  if(fmap) {
    fmap->mutex.lock();
    (fmap->refcount)--;
    ODINLOG(odinlog,normalDebug) << "fd/offset/refcount=" << fmap->fd << "/" << fmap->offset << "/" << fmap->refcount << STD_endl;
    if(!(fmap->refcount)) {
      fileunmap(fmap->fd, Array<T,N_rank>::data(), (LONGEST_INT)Array<T,N_rank>::size()*sizeof(T), fmap->offset);
      fmap->mutex.unlock();
      delete fmap;
      fmap=0;
    }
    if(fmap) fmap->mutex.unlock();
  }
}


template <typename T, int N_rank>
Data<T,N_rank>::~Data() {
  detach_fmap();
}


template <typename T, int N_rank>
void Data<T,N_rank>::reference(const Data<T,N_rank>& d) {
  Log<OdinData> odinlog("Data","reference");
  detach_fmap();
  fmap=d.fmap;
  if(fmap) {
    MutexLock lock(fmap->mutex);
    (fmap->refcount)++;
    ODINLOG(odinlog,normalDebug) << "fmap->refcount=" << fmap->refcount << STD_endl;
  }
  Array<T,N_rank>::reference(d);
}


template <typename T, int N_rank>
Data<T,N_rank>& Data<T,N_rank>::operator = (const tjarray<tjvector<T>,T>& a) {
  Log<OdinData> odinlog("Data","=");
  if( a.dim() == N_rank ) {
    ndim nn=a.get_extent();
    TinyVector<int,N_rank> tv;
    for(unsigned int i=0;i<a.dim();i++) tv(i)=nn[i];
    resize(tv);
    for(unsigned int i=0;i<a.total();i++) (*this)(create_index(i))=a[i];
  } else ODINLOG(odinlog,errorLog) << "dimension mismatch: this=" << N_rank<< ", tjarray=" << a.dim() << STD_endl;
  return *this;
}



template<typename T,int N_rank>
int Data<T,N_rank>::read(const STD_string& format,const STD_string& filename, LONGEST_INT offset) {
  Log<OdinData> odinlog("Data","read");
  if(format==TypeTraits::type2label((u8bit)0))  return read<u8bit> (filename);
  if(format==TypeTraits::type2label((s8bit)0))  return read<s8bit> (filename);
  if(format==TypeTraits::type2label((u16bit)0)) return read<u16bit>(filename);
  if(format==TypeTraits::type2label((s16bit)0)) return read<s16bit>(filename);
  if(format==TypeTraits::type2label((u32bit)0)) return read<u32bit>(filename);
  if(format==TypeTraits::type2label((s32bit)0)) return read<s32bit>(filename);
  if(format==TypeTraits::type2label((float)0))  return read<float> (filename);
  if(format==TypeTraits::type2label((double)0)) return read<double>(filename);
  ODINLOG(odinlog,errorLog) << "Unable to read file " << filename << " with data type " << format << STD_endl;
  return -1;
}



template<typename T,int N_rank>
template<typename T2>
int Data<T,N_rank>::read(const STD_string& filename, LONGEST_INT offset) {
  Log<OdinData> odinlog("Data","read");
  LONGEST_INT fsize=filesize(filename.c_str())-offset;
  LONGEST_INT nelements_file=fsize/sizeof(T2);
  LONGEST_INT length=product(Array<T,N_rank>::shape());
  ODINLOG(odinlog,normalDebug) << "fsize / nelements_file / length = " << fsize << " / " << nelements_file << " / " << length << STD_endl;

  if(!length) return 0;

  if(nelements_file<length) {
    ODINLOG(odinlog,errorLog) << "Size of file " << filename << " to small for reading" << STD_endl;
    return -1;
  }

  STD_string srctype=TypeTraits::type2label((T2)0);
  STD_string dsttype=TypeTraits::type2label((T)0);
  ODINLOG(odinlog,normalDebug) << "srctype/dsttype=" << srctype << "/" << dsttype << STD_endl;

  TinyVector<int,N_rank> fileshape(Array<T,N_rank>::shape());
  fileshape(N_rank-1) *= (Converter::get_elements((T)0)/Converter::get_elements((T2)0)); // Adjust extent of last dim for copmlex data
  Data<T2,N_rank> filedata(filename, true, fileshape, offset);
  filedata.convert_to(*this);

  return 0;
}



template<typename T,int N_rank>
int Data<T,N_rank>::write(const STD_string& format, const STD_string& filename, autoscaleOption scaleopt) const {
  Log<OdinData> odinlog("Data","write");
  if(format==TypeTraits::type2label((u8bit)0))  return write<u8bit> (filename,scaleopt);
  if(format==TypeTraits::type2label((s8bit)0))  return write<s8bit> (filename,scaleopt);
  if(format==TypeTraits::type2label((u16bit)0)) return write<u16bit>(filename,scaleopt);
  if(format==TypeTraits::type2label((s16bit)0)) return write<s16bit>(filename,scaleopt);
  if(format==TypeTraits::type2label((u32bit)0)) return write<u32bit>(filename,scaleopt);
  if(format==TypeTraits::type2label((s32bit)0)) return write<s32bit>(filename,scaleopt);
  if(format==TypeTraits::type2label((float)0))  return write<float> (filename,scaleopt);
  if(format==TypeTraits::type2label((double)0)) return write<double>(filename,scaleopt);
  ODINLOG(odinlog,errorLog) << "Unable to write file " << filename << " with data type " << format << STD_endl;
  return -1;
}


template<typename T,int N_rank>
int Data<T,N_rank>::write(const STD_string& filename, fopenMode mode) const {
  Log<OdinData> odinlog("Data","write");
  if(filename=="") return 0;

  FILE* file_ptr=FOPEN(filename.c_str(),modestring(mode));

  if(file_ptr==NULL) {
    ODINLOG(odinlog,errorLog) << "unable to create/open file >" << filename << "< - " << lasterr() << STD_endl;
    return -1;
  }

  Data<T,N_rank> data_copy(*this); // for contig memory

  LONGEST_INT nmemb=Array<T,N_rank>::numElements();
  LONGEST_INT count=fwrite(data_copy.c_array(),sizeof(T),nmemb,file_ptr);
  if(count!=nmemb) {
    ODINLOG(odinlog,errorLog) << "unable to fwrite to file >" << filename << "< - " << lasterr() << STD_endl;
    return -1;
  }

  if(file_ptr!=NULL) fclose(file_ptr);
  return 0;
}


template<typename T,int N_rank>
template<typename T2>
  int Data<T,N_rank>::write(const STD_string& filename, autoscaleOption scaleopt) const {
  Log<OdinData> odinlog("Data","write");

  rmfile(filename.c_str()); // remove old file to get new file size with mmap

  ODINLOG(odinlog,normalDebug) << "mem(pre): " << Profiler::get_memory_usage() << STD_endl;
  Data<T2,N_rank> converted_data; convert_to(converted_data,scaleopt);
  ODINLOG(odinlog,normalDebug) << "mem(cnv): " << Profiler::get_memory_usage() << STD_endl;
  Data<T2,N_rank> filedata(filename,false,converted_data.shape());
  ODINLOG(odinlog,normalDebug) << "mem(map): " << Profiler::get_memory_usage() << STD_endl;
  filedata=converted_data;

  return 0;
}



template<typename T,int N_rank>
int Data<T,N_rank>::write_asc_file(const STD_string& filename, const Array<T,N_rank>& pre, const Array<T,N_rank>& post) const {
  bool have_pre=false;
  bool have_post=false;

  Data<T,N_rank> pre_data(pre);
  Data<T,N_rank> post_data(post);
  int n=Array<T,N_rank>::numElements();

  if(pre_data.numElements()==n) have_pre=true;
  if(post_data.numElements()==n) have_post=true;

  STD_ofstream ofs(filename.c_str());
  if(ofs.bad()) return -1;

  T val;

  for(int i=0; i<n; i++) {
    if(have_pre) {
      val=pre_data(pre_data.create_index(i));
      ofs << val << " ";
    }
    val=(*this)(create_index(i));
    ofs << val;
    if(have_post) {
      val=post_data(post_data.create_index(i));
      ofs << " " << val;
    }
    ofs << "\n";
  }

  ofs.close();
  return 0;
}


template<typename T,int N_rank>
int Data<T,N_rank>::read_asc_file(const STD_string& filename) {
  STD_ifstream ifs(filename.c_str());
  if(ifs.bad()) return -1;

  STD_string valstr;
  for(int i=0; i<Array<T,N_rank>::numElements(); i++) {
    if(ifs.bad()) return -1;
    ifs >> valstr;
    TypeTraits::string2type(valstr,(*this)(create_index(i)));
  }

  ifs.close();
  return 0;
}



// Interface to FileIO functionality
int fileio_autowrite(const Data<float,4>& data, const STD_string& filename, const FileWriteOpts& opts, const Protocol* prot);
int fileio_autoread(Data<float,4>& data, const STD_string& filename, const FileReadOpts& opts, Protocol* prot);


template<typename T,int N_rank>
int Data<T,N_rank>::autowrite(const STD_string& filename, const FileWriteOpts& opts, const Protocol* prot) const {
  Data<float,4> filedata; convert_to(filedata);
  return fileio_autowrite(filedata, filename, opts, prot);
}

template<typename T,int N_rank>
int Data<T,N_rank>::autoread(const STD_string& filename, const FileReadOpts& opts, Protocol* prot) {
  Data<float,4> filedata;
  int result=fileio_autoread(filedata, filename, opts, prot);
  if(result>0) filedata.convert_to(*this);
  return result;
}



template<typename T,int N_rank>
Data<T,N_rank>::operator tjarray<tjvector<T>,T> () const {
   tjarray<tjvector<T>,T>  a;
   ndim nn(N_rank);
   for(unsigned int i=0; i<N_rank; i++) nn[i]=Array<T,N_rank>::extent(i);
   a.redim(nn);
   for(unsigned int i=0;i<a.total();i++) a[i]=(*this)(create_index(i));
   return a;
}




template<typename T,int N_rank>
template<typename T2, int N_rank2>
Data<T2,N_rank2>& Data<T,N_rank>::convert_to(Data<T2,N_rank2>& dst, autoscaleOption scaleopt) const {
  Log<OdinData> odinlog("Data","convert_to");

  TinyVector<int,N_rank2> newshape; newshape=1;
  for(int i=0; i<(N_rank-N_rank2+1); i++) { // collapse extra dims into first dim of new shape, including first dim
   int srcindex=i;
   if(srcindex>=0 && srcindex<N_rank) newshape(0)*=Array<T,N_rank>::extent(srcindex);
  }
  for(int i=0; i<(N_rank2-1); i++) { // Fill new shape with extents of original dim, except for the first dim
    int srcindex=N_rank-N_rank2+1+i;
    if(srcindex>=0 && srcindex<N_rank) newshape(1+i)=Array<T,N_rank>::extent(srcindex);
  }

  // modify last dimension to account for different element sizes
  newshape(N_rank2-1)=newshape(N_rank2-1)*Converter::get_elements((T)0)/Converter::get_elements((T2)0);

  dst.resize(newshape); // dst will be unique now

  Data<T,N_rank> src_copy(*this); // make read/write copy of this

  Converter::convert_array(src_copy.c_array(), dst.c_array(), src_copy.numElements(), dst.numElements(), scaleopt);

  return dst;
}


// specialization which uses reference if type/rank are same and no scaling is necessary to reduce memory usage
template<typename T,int N_rank>
Data<T,N_rank>& Data<T,N_rank>::convert_to(Data<T,N_rank>& dst, autoscaleOption scaleopt) const {
  Log<OdinData> odinlog("Data","convert_to");
  if(scaleopt==noscale || !std::numeric_limits<T>::is_integer) {
    ODINLOG(odinlog,normalDebug) << "Using reference" << STD_endl;
    dst.reference(*this);
  } else {
    Data<T,N_rank> src_copy(*this); // make read/write copy of this
    dst.resize(src_copy.shape()); // dst will be unique now
    Converter::convert_array(src_copy.c_array(), dst.c_array(), src_copy.numElements(), dst.numElements(), scaleopt);
  }
  return dst;
}


template<typename T,int N_rank>
TinyVector<int, N_rank> Data<T,N_rank>::create_index(unsigned long index) const {
  return index2extent<N_rank>(Array<T,N_rank>::shape(), index);
}


template<typename T,int N_rank>
unsigned long Data<T,N_rank>::create_linear_index(const TinyVector<int, N_rank>& indexvec) const {
  unsigned long totalIndex=0,subsize;
  TinyVector<int, N_rank> nn(Array<T,N_rank>::extent());
  for(unsigned long i=0;i<N_rank;i++) {
    nn(i)=1;
    subsize=product(nn);
    if(!subsize) subsize=1;
    totalIndex+=subsize*indexvec(i);
  }
  return totalIndex;
}


template<typename T,int N_rank>
void Data<T,N_rank>::shift(unsigned int shift_dim, int shift) {
  Log<OdinData> odinlog("Data","shift");

  if(!shift) return;

  if( shift_dim >= N_rank ){
    ODINLOG(odinlog,errorLog) << "shift dimension(" << shift_dim << ") >= rank of data (" << N_rank << ") !\n";
    return;
  }

  int shift_extent=Array<T,N_rank>::extent(shift_dim);
  int abs_shift=abs(shift);
  if(shift_extent < abs_shift){
    ODINLOG(odinlog,errorLog) << "extent(" << shift_extent << ") less than shift(" << abs_shift << ") !\n";
    return;
  }

  Data<T,N_rank> data_copy(Array<T,N_rank>::copy());

  TinyVector<int,N_rank> index;
  for(int i=0; i<Array<T,N_rank>::numElements(); i++) {
    index=create_index(i);
    T val=data_copy(index);
    int shiftindex=index(shift_dim)+shift;
    if(shiftindex>=shift_extent) shiftindex-=shift_extent;
    if(shiftindex<0)             shiftindex+=shift_extent;
    index(shift_dim)=shiftindex;
    (*this)(index)=val;
  }

}




template<typename T,int N_rank>
void Data<T,N_rank>::congrid(const TinyVector<int,N_rank>& newshape, const TinyVector<float,N_rank>* subpixel_shift, bool left_to_right) {
  Log<OdinData> odinlog("Data","congrid");
  for(int i=0;i<N_rank;i++) {
    ODINLOG(odinlog,normalDebug) << "oldshape/newshape(" << i << ")=" << Array<T,N_rank>::shape()(i) << "/" << newshape(i) << STD_endl;
    if(subpixel_shift) ODINLOG(odinlog,normalDebug) << "subpixel_shift(" << i << ")=" << (*subpixel_shift)(i) << STD_endl;
  }

  for(int irank=0; irank<N_rank; irank++) {
    int dim=irank;
    if(!left_to_right) dim=N_rank-1-irank;
    float shift=0.0;
    if(subpixel_shift) shift=(*subpixel_shift)(dim);
    interpolate1dim(dim,newshape(dim),shift);
  }
}



template<typename T,int N_rank>
void Data<T,N_rank>::interpolate1dim(unsigned int dim, int newsize, float subpixel_shift) {
  Log<OdinData> odinlog("Data","interpolate1dim");
  ODINLOG(odinlog,normalDebug) << "dim/newsize=" << dim << "/" << newsize << STD_endl;

  if(newsize==Array<T,N_rank>::shape()(dim) && subpixel_shift==0.0) {
    ODINLOG(odinlog,normalDebug) << "interpolation not required" << STD_endl;
    return;
  }

  if(dim>=N_rank) {
    ODINLOG(odinlog,errorLog) << "dim is larger than N_rank" << STD_endl;
    return;
  }
  if(newsize<0) {
    ODINLOG(odinlog,errorLog) << "newsize is negative" << STD_endl;
    return;
  }

  Array<T,N_rank> olddata(*this);
  olddata.makeUnique();

  TinyVector<int,N_rank> oldshape(olddata.shape());
  int oldsize=oldshape(dim);

  TinyVector<int,N_rank> newshape(oldshape);
  newshape(dim)=newsize;
  resize(newshape);

  // This holds the the shape of the subspace which is orthogonal to direction 'dim'
  TinyVector<int,N_rank> ortho_shape(oldshape);
  ortho_shape(dim)=1;

  // The total number of elements in the orthogonal subspace
  unsigned long n_ortho=product(ortho_shape);

  ODINLOG(odinlog,normalDebug) << "n_ortho=" << n_ortho << STD_endl;

  TinyVector<int,N_rank> indexvec;
  T* oldoneline=new T[oldsize];

  for(unsigned long iortho=0; iortho<n_ortho; iortho++) {

    indexvec=index2extent<N_rank>(ortho_shape,iortho);

    for(int j=0; j<oldsize; j++) {
      indexvec(dim)=j;
      oldoneline[j]=olddata(indexvec);
    }

    T* newoneline=interpolate1D(oldoneline,oldsize,newsize,subpixel_shift);

    for(int j=0; j<newsize; j++) {
      indexvec(dim)=j;
      (*this)(indexvec)=newoneline[j];
    }

    delete[] newoneline;
  }

  delete[] oldoneline;
}



template<typename T,int N_rank>
T* Data<T,N_rank>::c_array() {
  Log<OdinData> odinlog("Data","c_array");

  bool need_copying=false;

  // check storage order
  TinyVector<int,N_rank>ord=Array<T,N_rank>::ordering();
  for(int i=0; i<N_rank-1; i++) {
    if(ord(i)<ord(i+1)) need_copying=true;
  }

  // check storage direction
  for(int i=0; i<N_rank; i++) {
    if(!Array<T,N_rank>::isRankStoredAscending(i)) need_copying=true;
  }

  // check for slicing
  if(!Array<T,N_rank>::isStorageContiguous()) need_copying=true;

  if(need_copying) {
    ODINLOG(odinlog,normalDebug) << "need_copying" << STD_endl;
    Data<T,N_rank> tmp(Array<T,N_rank>::shape());
    tmp=(*this); // reference(copy()) would not be enough since it would preserve ordering
    reference(tmp);
  }

  return Array<T,N_rank>::data();
}

#endif

---