source: Sophya/trunk/SophyaLib/TArray/tarray.h@ 1517

// 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 "basarr.h"
#include "ndatablock.h"
#include <complex>
#include "utilarr.h"


namespace SOPHYA {

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

//   --------------------------- classe template Array -----------------------
//  ( See BaseArray class for data organisation in  memory and related methods )

//! Class for template arrays
template <class T>
class TArray : public BaseArray {
public:
  // Creation / destruction 
  TArray();
  TArray(int_4 ndim, const sa_size_t * siz, sa_size_t step =1);
  TArray(sa_size_t nx, sa_size_t ny=0, sa_size_t nz=0, sa_size_t nt=0, sa_size_t nu=0);
  TArray(int_4 ndim, const sa_size_t * siz, NDataBlock<T> & db, bool share=false, sa_size_t step=1, sa_size_t offset=0);
  TArray(int_4 ndim, const sa_size_t * siz, T* values, sa_size_t step=1, sa_size_t offset=0, Bridge* br=NULL);
  TArray(const TArray<T>& a);
  TArray(const TArray<T>& a, bool share);
  TArray(const BaseArray& a);

  virtual ~TArray();

  // A = B
  //! = operator between TArray
  /*! \warning Datas are copied (cloned) from \b a.
      \sa Set \sa NDataBlock::operator=(const NDataBlock<T>&) */
  inline  TArray<T>& operator = (const TArray<T>& a) { return Set(a); }
  virtual TArray<T>& Set(const TArray<T>& a);

  //! = operator between  TArray 's with different types - Elements are converted.
  inline  TArray<T>& operator = (const BaseArray& a) { return SetBA(a); }
  virtual TArray<T>& SetBA(const BaseArray& a);

  // Gestion taille/Remplissage
  virtual void Clone(const TArray<T>& a); 
  // 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);

  void ReSize(int_4 ndim, sa_size_t * siz, sa_size_t step=1);
  void ReSize(const BaseArray& a);
  //! a synonym (alias) for method ReSize(int_4, ...)
  inline void SetSize(int_4 ndim, sa_size_t * siz, sa_size_t step=1)
                { ReSize(ndim, siz, step); }
  //! a synonym (alias) for method ReSize(const BaseArray&)
  inline void SetSize(const BaseArray& a)
                { ReSize(a); }
  void Realloc(int_4 ndim, sa_size_t * siz, sa_size_t step=1, bool force=false);

  // Compacts size=1 array dimensions
  virtual TArray<T>& CompactAllDimensions();
  virtual TArray<T>& CompactTrailingDimensions();

  // Packing array elements in memory
  virtual TArray<T> PackElements(bool force=false) const ;

  // SubArrays - $CHECK$ Reza 03/2000 je ne sais pas s'il faut declarer ca const ??
  TArray<T> SubArray(Range rx, Range ry, Range rz, Range rt, Range ru) const ;
#if !defined(__GNUG__)
  //! () operator for Sub arrays extraction
  /*! \sa SubArray */
  inline TArray<T> operator () (Range rx, Range ry, Range rz, Range rt=0, Range ru=0) const
                   { return  SubArray(rx, ry, rz, rt, ru); }
#else
  // g++ (2.95) se melange les pinceaux avec les arguments par defaut !
  //! () operator for Sub arrays extraction
  /*! \sa SubArray */
  inline TArray<T> operator () (Range rx, Range ry, Range rz) const
                   { return  SubArray(rx, ry, rz, Range(0), Range(0)); }
  inline TArray<T> operator () (Range rx, Range ry, Range rz, Range rt) const
                   { return  SubArray(rx, ry, rz, rt, Range(0)); }
  inline TArray<T> operator () (Range rx, Range ry, Range rz, Range rt, Range ru) const
                   { return  SubArray(rx, ry, rz, rt, ru); }
#endif

  // ---- Access to data
  // Definition of virtual element acces method inherited from BaseArray class 
  virtual MuTyV & ValueAtPosition(sa_size_t ip) const;

  // Data Access: operator overloaded inline acces methods 
  inline T const& operator()(sa_size_t ix, sa_size_t iy, sa_size_t iz) const ;
  inline T&       operator()(sa_size_t ix, sa_size_t iy, sa_size_t iz);
  inline T const& operator()(sa_size_t ix, sa_size_t iy, sa_size_t iz, sa_size_t it, sa_size_t iu=0) const ;
  inline T&       operator()(sa_size_t ix, sa_size_t iy, sa_size_t iz, sa_size_t it, sa_size_t iu=0);
  inline T const& operator[](sa_size_t ip) const ;
  inline T&       operator[](sa_size_t ip);

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

  //! Return pointer to first element adress
  inline        T* Data()       {return mNDBlock.Begin()+offset_;}
  //! Return pointer to first element adress
  inline const  T* Data() const {return mNDBlock.Begin()+offset_;}
  //! Return reference to datablock NDataBlock
  inline        NDataBlock<T>& DataBlock()       {return mNDBlock;}
  //! Return reference to datablock NDataBlock
  inline const  NDataBlock<T>& DataBlock() const {return mNDBlock;}

  // Temporaire?
  //! Are the array temporay ?
  inline bool   IsTemp(void) const {return mNDBlock.IsTemp();}
  //! Set the array as temporay
  inline void   SetTemp(bool temp=false) const {mNDBlock.SetTemp(temp);}

// Operations diverses  = , +=, ...
// Conversion en type T, if Size() == 1
  inline T toScalar();
// Met les elements a une suite de valeurs
  virtual TArray<T>&  SetSeq(Sequence const & seq);
  //! Fill TArray with Sequence \b seq
  inline  TArray<T>&  operator = (Sequence const & seq)    { return SetSeq(seq); }
// A = x (tous les elements a x)
  virtual TArray<T>&  SetT(T x);
  //! Fill TArray with all elements equal to \b x
  inline  TArray<T>&  operator = (T x)             { return SetT(x); }
// A += -= *= /= x (ajoute, soustrait, ... x a tous les elements)
  virtual TArray<T>&  Add(T x);
  //! Add \b x to all elements
  inline  TArray<T>&  operator += (T x)            { return Add(x); }
  virtual TArray<T>&  Sub(T x, bool fginv=false);
  //! Substract \b x to all elements
  inline  TArray<T>&  operator -= (T x)            { return Sub(x); }
  virtual TArray<T>&  Mul(T x);
  //! Multiply all elements by \b x
  inline  TArray<T>&  operator *= (T x)            { return Mul(x); }
  virtual TArray<T>&  Div(T x, bool fginv=false);
  //! Divide all elements by \b x
  inline  TArray<T>&  operator /= (T x)            { return Div(x); }

// applique le signe moins a tous les elements
  virtual TArray<T>&  NegateElt();
  
// A += -=  (ajoute, soustrait element par element les deux tableaux )
  virtual TArray<T>&  AddElt(const TArray<T>& a);
  //! Operator TArray += TArray
  inline  TArray<T>&  operator += (const TArray<T>& a)  { return AddElt(a); }
  virtual TArray<T>&  SubElt(const TArray<T>& a, bool fginv=false);
  //! Operator TArray -= TArray
  inline  TArray<T>&  operator -= (const TArray<T>& a)  { return SubElt(a); }
// Multiplication, division element par element les deux tableaux 
  virtual TArray<T>&  MulElt(const TArray<T>& a);
  virtual TArray<T>&  DivElt(const TArray<T>& a, bool fginv=false, bool divzero=false);
// Recopie des valeurs, element par element
  virtual TArray<T>&  CopyElt(const TArray<T>& a);
// Recopie des valeurs avec conversion prealable, element par element
  virtual TArray<T>&  ConvertAndCopyElt(const BaseArray& a);

// Somme et produit des elements
  virtual T Sum() const ;
  virtual T Product() const ;
// Somme du carre des elements
  virtual T SumX2() const;
// Valeur min et max des elements
  virtual void MinMax(T& min, T& max) const ;

// Impression, I/O, ...
  virtual string InfoString() const;   
  virtual void   Print(ostream& os, int_4 maxprt=-1, bool si=false) const ;

// Lecture,Ecriture sur fichier ASCII
  virtual void   ReadASCII(istream& is);
  virtual void   WriteASCII(ostream& os) const;

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

protected:

  NDataBlock<T> mNDBlock; //!< Block for datas
  mutable MuTyV my_mtv;   //!< for use by ValueAtPosition()
};

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

//   Surcharge d'operateur >> 
//! Decodes the ASCII input stream \b is , filling TArray \b a elements
template <class T>
inline istream& operator >> (istream& is, TArray<T>& a)
                            { a.ReadASCII(os);    return(is);    }


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

/*! \ingroup TArray \fn operator+(const TArray<T>&,T)
  \brief Operator TArray = TArray + constant */
template <class T> inline TArray<T> operator + (const TArray<T>& a, T b)
    {TArray<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Add(b); return result;}

/*! \ingroup TArray \fn operator+(T,const TArray<T>&)
  \brief Operator TArray = constant + TArray */
template <class T> inline TArray<T> operator + (T b,const TArray<T>& a)
    {TArray<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Add(b); return result;}

/*! \ingroup TArray \fn operator-(const TArray<T>&,T)
  \brief Operator TArray = TArray - constant */
template <class T> inline TArray<T> operator - (const TArray<T>& a, T b)
    {TArray<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Sub(b); return result;}

/*! \ingroup TArray \fn operator-(T,const TArray<T>&)
  \brief Operator TArray = constant - TArray */
template <class T> inline TArray<T> operator - (T b,const TArray<T>& a)
    {TArray<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Sub(b,true); return result;}

/*! \ingroup TArray \fn operator*(const TArray<T>&,T)
  \brief Operator TArray = TArray * constant */
template <class T> inline TArray<T> operator * (const TArray<T>& a, T b)
    {TArray<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Mul(b); return result;}

/*! \ingroup TArray \fn operator*(T,const TArray<T>&)
  \brief Operator TArray = constant * TArray */
template <class T> inline TArray<T> operator * (T b,const TArray<T>& a)
    {TArray<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Mul(b); return result;}

/*! \ingroup TArray \fn operator/(const TArray<T>&,T)
  \brief Operator TArray = TArray / constant */
template <class T> inline TArray<T> operator / (const TArray<T>& a, T b)
    {TArray<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Div(b); return result;}

/*! \ingroup TArray \fn operator/(T,const TArray<T>&)
  \brief Operator TArray = constant / TArray  */
template <class T> inline TArray<T> operator / (T b, const TArray<T>& a)
    {TArray<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Div(b, true); return result;}

////////////////////////////////////////////////////////////////
// Surcharge d'operateurs B = -A

/*! \ingroup TArray \fn operator - (const TArray<T>&)
  \brief Operator - Returns an array with elements equal to the opposite of
  the original array elements.  */
template <class T> inline TArray<T> operator - (const TArray<T>& a)
    {TArray<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.NegateElt(); return result;}

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

/*! \ingroup TArray \fn operator+(const TArray<T>&,const TArray<T>&)
  \brief Operator TArray = TArray + TArray */
template <class T>
inline TArray<T> operator + (const TArray<T>& a,const TArray<T>& b)
    { TArray<T> result; result.SetTemp(true); 
    if (b.IsTemp())  { result.Share(b); result.AddElt(a); }
    else { result.CloneOrShare(a); result.AddElt(b); }
    return result; }

/*! \ingroup TArray \fn operator-(const TArray<T>&,const TArray<T>&)
  \brief Operator TArray = TArray - TArray */
template <class T>
inline TArray<T> operator - (const TArray<T>& a,const TArray<T>& b)
    { TArray<T> result; result.SetTemp(true); 
    if (b.IsTemp())  { result.Share(b); result.SubElt(a, true); }
    else { result.CloneOrShare(a); result.SubElt(b); }
    return result; }


// --------------------------------------------------
//        inline element acces methods
// --------------------------------------------------

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

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

//! Return element (ix,iy,iz,it,iu) value with Check of indexes bound first
template <class T>
inline T const& TArray<T>::ElemCheckBound(sa_size_t ix, sa_size_t iy, sa_size_t iz, sa_size_t it, sa_size_t iu) const
{
  CheckBound(ix, iy, iz, it, iu, 4);
  return(Elem(ix, iy, iz, it, iu));
}

//! Return element (ix,iy,iz,it,iu) value with Check of indexes bound first
template <class T>
inline T & TArray<T>::ElemCheckBound(sa_size_t ix, sa_size_t iy, sa_size_t iz, sa_size_t it, sa_size_t iu)
{
  CheckBound(ix, iy, iz, it, iu, 4);
  return(Elem(ix, iy, iz, it, iu));
}

//! Return element (ix,iy,iz) value
template <class T>
inline T const& TArray<T>::operator()(sa_size_t ix, sa_size_t iy, sa_size_t 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]) );
}

//! Return element (ix,iy,iz) value
template <class T>
inline T & TArray<T>::operator()(sa_size_t ix, sa_size_t iy, sa_size_t 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]) );
}

//! Operator () : return element (ix,iy,iz,it,iu) value
template <class T>
inline T const& TArray<T>::operator()(sa_size_t ix, sa_size_t iy, sa_size_t iz, sa_size_t it, sa_size_t 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]) );
}

//! Operator () : return element (ix,iy,iz,it,iu) value
template <class T>
inline T & TArray<T>::operator()(sa_size_t ix, sa_size_t iy, sa_size_t iz, sa_size_t it, sa_size_t 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]) );
}


//! Operator [] : return element at positon ip
template <class T>
inline T const& TArray<T>::operator[](sa_size_t ip) const
{
#ifdef SO_BOUNDCHECKING
  if (ip >= totsize_)  throw( ParmError("TArray<T>::operator[] Out-of-bound Error") );
#endif
return *(mNDBlock.Begin()+Offset(ip));
}

//! Operator [] : return element at positon ip
template <class T>
inline T & TArray<T>::operator[](sa_size_t ip) 
{
#ifdef SO_BOUNDCHECKING
  if (ip >= totsize_)  throw( ParmError("TArray<T>::operator[] Out-of-bound Error") );
#endif
return *(mNDBlock.Begin()+Offset(ip));
}


//! Converts to a scalar (value of first element) if the array size is equal to 1 
template <class T>
inline T TArray<T>::toScalar()
{ 
  if (Size() != 1) throw(SzMismatchError("TArray<T>::operator T() Size() != 1")) ;
  return ( (*this)[0] );
}

// Typedef pour simplifier
/*! \ingroup TArray
  \typedef Array
  \brief To simplified TArray<r_8> writing
*/
typedef TArray<r_8> Array;

} // Fin du namespace

#endif