// This may look like C code, but it is really -*- C++ -*-
//      template array class for numerical types
//                     R. Ansari, C.Magneville   03/2000

#ifndef TArray_SEEN
#define TArray_SEEN

#include "machdefs.h"
#include <math.h>
#include <iostream.h>
#include "ndatablock.h"
#include <complex>
#include "anydataobj.h"
#include "ndatablock.h"
#include "dvlist.h"


//  Maximum number of dimensions for array
#define TARRAY_MAXNDIMS  5

namespace SOPHYA {

/* Quelques utilitaires */
typedef double (* Arr_DoubleFunctionOfX) (double x);
typedef float  (* Arr_FloatFunctionOfX)  (float x);

class Sequence { 
public:
  Sequence (double start=0., double step=1.);
  inline double & Start() { return start_; }
  inline double & Step() { return step_; }
  inline double operator () (uint_4 k) { return(start_+(double)k*step_); }
protected:
  double start_, step_;
};

// Forward declaration
template <class T> class FIO_TArray;

//   ------------ classe template Array ----------- 
template <class T>
class TArray : public AnyDataObj {
public:
  // Creation / destruction 
  TArray();
  TArray(uint_4 ndim, uint_4 * siz, uint_4 step =1);
  TArray(uint_4 nx, uint_4 ny=0, uint_4 nz=0);
  TArray(uint_4 ndim, uint_4 * siz, NDataBlock<T> & db, bool share=false, uint_4 step=1, uint_8 offset=0);
  TArray(uint_4 ndim, uint_4 * siz, T* values, uint_4 step=1, uint_8 offset=0, Bridge* br=NULL);
  TArray(const TArray<T>& a);
  TArray(const TArray<T>& a, bool share);

  virtual ~TArray();

  // A = B : partage les donnees si "a" est temporaire, clone sinon.
  virtual TArray<T>& operator = (const TArray<T>& a);

  // Extraction de sous-tableau
  virtual TArray<T> SubArray(uint_4 ndim, uint_4 * siz, uint_4 * pos);

  // Gestion taille/Remplissage
  virtual void Clone(const TArray<T>& a); 
  virtual void ReSize(uint_4 ndim, uint_4 * siz, uint_4 step=1);
  virtual void Realloc(uint_4 ndim, uint_4 * siz, uint_4 step=1, bool force=false);
  // Returns true if ndim and sizes are equal
  virtual bool CompareSizes(const TArray<T>& a);
 
  // Array dimensions
  inline uint_4 NbDimensions() const { return( ndim_ ); }

  inline uint_8 Size() const { return(totsize_); }
  inline uint_4 SizeX() const { return(size_[0]); }
  inline uint_4 SizeY() const { return(size_[1]); }
  inline uint_4 SizeZ() const { return(size_[2]); }
  inline uint_4 Size(int ka) const { return(size_[CheckDI(ka,1)]); }

  // memory organisation - packing information
  inline bool   IsPacked() const { return(minstep_ == 1); }
  inline bool   IsPackedX() const { return(step_[0] == 1); }
  inline bool   IsPackedY() const { return(step_[1] == 1); }
  inline bool   IsPackedZ() const { return(step_[2] == 1); }
  inline bool   IsPacked(int ka) const { return(step_[CheckDI(ka,2)] == 1); }

  inline uint_4 Step() const  { return(minstep_); }
  inline uint_4 StepX() const { return(step_[0]); }
  inline uint_4 StepZ() const { return(step_[1]); }
  inline uint_4 StepY() const { return(step_[2]); }
  inline uint_4 Step(int ka) const { return(step_[CheckDI(ka,3)]); }

  inline uint_8 Offset() const { return(offset_); }


  // Temporaire?
  inline bool   IsTemp(void) const {return mNDBlock.IsTemp();}
  inline void   SetTemp(bool temp=false) const {mNDBlock.SetTemp(temp);}

  // Acces to data
  inline T const& operator()(uint_4 ix, uint_4 iy, uint_4 iz) const ;
  inline T&       operator()(uint_4 ix, uint_4 iy, uint_4 iz);
  inline T const& operator()(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu=0) const ;
  inline T&       operator()(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu=0);
  inline T const& operator[](uint_4 ip) const ;
  inline T&       operator[](uint_4 ip);

  inline T const& Elem(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it=0, uint_4 iu=0) const ;
  inline T&       Elem(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it=0, uint_4 iu=0);
  inline T const& ElemCheckBound(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it=0, uint_4 iu=0) const ;
  inline T&       ElemCheckBound(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it=0, uint_4 iu=0);

  inline        T* Data()       {return mNDBlock.Begin()+offset_;}
  inline const  T* Data() const {return mNDBlock.Begin()+offset_;}
  inline        NDataBlock<T>& DataBlock()       {return mNDBlock;}
  inline const  NDataBlock<T>& DataBlock() const {return mNDBlock;}

// Operations diverses  = , +=, ...
// Conversion en type T, if Size() == 1
  inline operator T();
// Met les elements a une suite de valeurs
  virtual TArray<T>&  operator = (Sequence & seq);
// A = x (tous les elements a x)
  virtual TArray<T>&  operator = (T x);
// A += -= *= /= x (ajoute, soustrait, ... x a tous les elements)
  virtual TArray<T>&  operator += (T x);
  virtual TArray<T>&  operator -= (T x);
  virtual TArray<T>&  operator *= (T x);
  virtual TArray<T>&  operator /= (T x);
// A += -=  (ajoute, soustrait element par element les deux tableaux )
  virtual TArray<T>&  operator += (const TArray<T>& a);
  virtual TArray<T>&  operator -= (const TArray<T>& a);
// Multilication, division element par element les deux tableaux 
  virtual TArray<T>&  MultElt(const TArray<T>& a);
  virtual TArray<T>&  DivElt(const TArray<T>& a);

// Application d'une fonction
  virtual TArray<T>&  ApplyFunction(Arr_DoubleFunctionOfX f);
  virtual TArray<T>&  ApplyFunction(Arr_FloatFunctionOfX f);

// Impression, I/O, ...
  static void   SetMaxPrint(uint_4 nprt=50);   
  void          Show(ostream& os, bool si=false) const;
  inline void   Show() const { Show(cout); }
  void          Print(ostream& os, int_4 maxprt=-1, bool si=false) const ;

//  Objet DVList info 
  DVList&       Info();

//  Pour la gestion de persistance
  friend class  FIO_TArray<T>;

protected:
  // Verifie la compatibilite de l'index de dimension
  inline int CheckDI(int ka, int msg) const ;
  // Verifie la compatibilite des bornes d'index
  inline void CheckBound(int ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu, int msg) const ;
  // Changing Sizes/NDim ... return true if OK
  bool UpdateSizes(uint_4 ndim, const uint_4 * siz, uint_4 step, uint_8 offset, string & exmsg);
  bool UpdateSizes(uint_4 ndim, const uint_4 * siz, const uint_4 * step, uint_8 offset, string & exmsg);
  bool UpdateSizes(const TArray<T>& a, string & exmsg);
  static uint_8 ComputeTotalSize(uint_4 ndim, uint_4 * siz, uint_4 step, uint_8 offset) ;
  // partage les donnees si "a" temporaire, clone sinon.
  void CloneOrShare(const TArray<T>& a);
  // Share: partage les donnees de "a"
  void Share(const TArray<T>& a);

  uint_4 ndim_;                   // nb of dimensions
  uint_4 size_[TARRAY_MAXNDIMS];     // array size in each dimension
  uint_8 totsize_;                   // Total number of elements
  uint_4 step_[TARRAY_MAXNDIMS];     // two consecutive elements distance in a given dimension 
  uint_4 minstep_;                   // minimal step (in any axes)
  uint_8 offset_;              // global offset -> position of elem[0] in DataBlock
  NDataBlock<T> mNDBlock;      // Le bloc des donnees
  uint_4 marowi_, macoli_;     // For matrices, Row index and column index in dimensions

  DVList* mInfo;               // Infos (variables) attachees au tableau

  static char * ck_op_msg_[6];  // Operation messages for CheckDI() CheckBound()
  static uint_4 max_nprt_;      // Nb maxi d'elements imprimes
};

////////////////////////////////////////////////////////////////
//   Surcharge d'operateur << 
template <class T>
inline ostream& operator << (ostream& os, const TArray<T>& a)
                            { a.Print(os);    return(os);    }

////////////////////////////////////////////////////////////////
// Surcharge d'operateurs A (+,-,*,/) (T) x

template <class T> inline TArray<T> operator + (const TArray<T>& a, T b)
    {TArray<T> result(a); result.SetTemp(true); result += b; return result;}

template <class T> inline TArray<T> operator + (T b,const TArray<T>& a)
    {TArray<T> result(a); result.SetTemp(true); result += b; return result;}

template <class T> inline TArray<T> operator - (const TArray<T>& a, T b)
    {TArray<T> result(a); result.SetTemp(true); result -= b; return result;}

template <class T> inline TArray<T> operator - (T b,const TArray<T>& a)
    {TArray<T> result(a); result.SetTemp(true);
     result.DataBlock() = b-result.DataBlock(); return result;}

template <class T> inline TArray<T> operator * (const TArray<T>& a, T b)
    {TArray<T> result(a); result.SetTemp(true); result *= b; return result;}

template <class T> inline TArray<T> operator * (T b,const TArray<T>& a)
    {TArray<T> result(a); result.SetTemp(true); result *= b; return result;}

template <class T> inline TArray<T> operator / (const TArray<T>& a, T b)
    {TArray<T> result(a); result.SetTemp(true); result /= b; return result;}

////////////////////////////////////////////////////////////////
// Surcharge d'operateurs C = A (+,-) B

template <class T>
inline TArray<T> operator + (const TArray<T>& a,const TArray<T>& b)
    {TArray<T> result(a); result.SetTemp(true); result += b; return result;}

template <class T>
inline TArray<T> operator - (const TArray<T>& a,const TArray<T>& b)
    {TArray<T> result(a); result.SetTemp(true); result += b; return result;}

////////////////////////////////////////////////
//  Calcul de fonction 
template <class T>
inline TArray<T> sin(const TArray<T>& a)
    {TArray<T> result(a); result.SetTemp(true); result.ApplyFunction(sin); return result;}

template <class T>
inline TArray<T> cos(const TArray<T>& a)
    {TArray<T> result(a); result.SetTemp(true); result.ApplyFunction(cos); return result;}

// --------------------------------------------------
//        Methodes inline de verification
// --------------------------------------------------
template <class T>
inline int TArray<T>::CheckDI(int ka, int msg)  const
{
  if ( (ka < 0) || (ka >= ndim_) ) {
    string txt = "TArray<T>::CheckDimensionIndex/Error ";   txt += ck_op_msg_[msg];
    throw(ParmError(txt));
  }
  return(ka);
}

template <class T>
inline void TArray<T>::CheckBound(int ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu, int msg)  const
{
  if ( (ix >= size_[0]) || (iy >= size_[1]) || (iz > size_[2]) ||
       (it >= size_[3]) || (iu >= size_[4]) ) {
    string txt = "TArray<T>::CheckArrayBound/Error ";   txt += ck_op_msg_[msg];
    throw(ParmError(txt));
  }
  return;
}

// --------------------------------------------------
//        inline element acces methods
// --------------------------------------------------
template <class T>
inline T const& TArray<T>::Elem(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu) const
{
  return ( *( mNDBlock.Begin()+ offset_+ 
	      ix*step_[0] + iy*step_[1] + iz*step_[2] +
	      it*step_[3] + iu*step_[4]) );
}

template <class T>
inline T & TArray<T>::Elem(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu)
{
  return ( *( mNDBlock.Begin()+ offset_+ 
	      ix*step_[0] + iy*step_[1] + iz*step_[2] +
	      it*step_[3] + iu*step_[4]) );
}

template <class T>
inline T const& TArray<T>::ElemCheckBound(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu) const
{
  CheckBound(ix, iy, iz, it, iu, 4);
  Elem(ix, iy, iz, it, iu);
}

template <class T>
inline T & TArray<T>::ElemCheckBound(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu)
{
  CheckBound(ix, iy, iz, it, iu, 4);
  Elem(ix, iy, iz, it, iu);
}

template <class T>
inline T const& TArray<T>::operator()(uint_4 ix, uint_4 iy, uint_4 iz) const
{
#ifdef SO_BOUNDCHECKING
  CheckBound(ix, iy, iz, 0, 0, 4);
#endif
  return ( *( mNDBlock.Begin()+ offset_+ 
	      ix*step_[0] + iy*step_[1] + iz*step_[2]) );
}

template <class T>
inline T & TArray<T>::operator()(uint_4 ix, uint_4 iy, uint_4 iz) 
{
#ifdef SO_BOUNDCHECKING
  CheckBound(ix, iy, iz, 0, 0, 4);
#endif
  return ( *( mNDBlock.Begin()+ offset_+ 
	      ix*step_[0] + iy*step_[1] + iz*step_[2]) );
}

template <class T>
inline T const& TArray<T>::operator()(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu) const
{
#ifdef SO_BOUNDCHECKING
  CheckBound(ix, iy, iz, it, iu, 4);
#endif
  return ( *( mNDBlock.Begin()+ offset_+ 
	      ix*step_[0] + iy*step_[1] + iz*step_[2] +
	      it*step_[3] + iu*step_[4]) );
}

template <class T>
inline T & TArray<T>::operator()(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu)
{
#ifdef SO_BOUNDCHECKING
  CheckBound(ix, iy, iz, it, iu, 4);
#endif
  return ( *( mNDBlock.Begin()+ offset_+ 
	      ix*step_[0] + iy*step_[1] + iz*step_[2] +
	      it*step_[3] + iu*step_[4]) );
}

template <class T>
inline T const& TArray<T>::operator[](uint_4 ip) const
{
#ifdef SO_BOUNDCHECKING
  if (ip >= totsize_)  throw( ParmError("TArray<T>::operator[] Out-of-bound Error") );
#endif
return *(mNDBlock.Begin()+offset_+ip*minstep_);
}

template <class T>
inline T & TArray<T>::operator[](uint_4 ip) 
{
#ifdef SO_BOUNDCHECKING
  if (ip >= totsize_)  throw( ParmError("TArray<T>::operator[] Out-of-bound Error") );
#endif
return *(mNDBlock.Begin()+offset_+ip*minstep_);
}

template <class T>
inline TArray<T>::operator T()
{ 
  if (Size() != 1) throw(SzMismatchError("TArray<T>::operator T() Size() != 1")) ;
  return ( (*this)[0] );
}

} // Fin du namespace

#endif