#include "sopnamsp.h"
#include "machdefs.h"

#include <math.h>
#include <iostream>
#include <fstream>

#include "tarrinit.h"
#include "array.h"
#include "timing.h"

/*  --------------------------------------------------------- */
/*  ------- Programme de test des fonctions de base  -------  */
/*  ----    sur les TArray/TMatrix de SOPHYA    ----          */
/*  2000-2006      -   C. Magneville,  R. Ansari              */
/*  arrt -> Help/Usage                                        */
/*  Pour effectuer les tests avec RC                          */
/*  arrt  check       -> Rc=0 OK                              */
/*  Pour verification avec les impressions                    */
/*  arrt a / w / m / ascio                                    */
/*  --------------------------------------------------------- */

static void tstarr(bool tio);
static void tstasciiio();
static void tstioarr(TArray<int_4> & ia, TArray<int_4> & ib, TArray<int_4> & ic);
static void tstmtx(int n);
void _ZZ_TestTMatrixRC_YY_(TMatrix<r_8> & m);
// Methode de verification des fonctionalites - Fev 2006 
static int checktarr();

static int prtlev = 1;

int main(int narg, char* arg[])
{

  SophyaInit();
  InitTim();   // Initializing the CPU timer


  if (narg < 2) {
    cout << " arrt - TArray<T> Test - Usage arrt a/w/m/z [size=5] [prtlev=1] [nprt=50]\n" 
	 << " a: Simple Array Test  w: a+ PPersist Test (FIO_TArray<T>) \n" 
	 << " m: Matrix and Vector Test \n" 
	 << " check: TArray Test/Check with Rc  \n" 
	 << " ascio: Test of ascii I/O \n" 
	 << " z: Appel_ZZ_TestTMatrixRC_YY_ " << endl;
    exit(0);
  }

  string opt = arg[1];
  int n = 5;
  if (narg > 2) n = atoi(arg[2]);
  uint_4 nprt = 50;
  if (narg > 3) prtlev = atoi(arg[3]);
  if (narg > 4) nprt = atoi(arg[4]);

  BaseArray::SetMaxPrint(nprt, prtlev);

  int rc = 0;
  try {
    if (opt == "w")  tstarr(true);
    else if (opt == "a") tstarr(false);
    else if (opt == "z") {
      if (n<5) n = 5;
      TMatrix<r_8> tza(n,n+3);
      tza = RegularSequence(2.,1.5);
      _ZZ_TestTMatrixRC_YY_(tza);
      TMatrix<r_8> tza2 = tza(Range(1,3), Range(2,5));
      _ZZ_TestTMatrixRC_YY_(tza2);
    }
    else if (opt == "m") tstmtx(n);
    else if (opt == "ascio") tstasciiio();
    else if (opt == "check") rc = checktarr();
    else {
      cout << " arrt/Error : unknown option ! " << endl;
    }
  }
  catch (PThrowable exc) {
    cerr << " catched Exception " << exc.Msg() << endl;
    rc = 97;
  }  
  catch (...) {
    cerr << " catched unknown (...) exception " << endl; 
    rc  = 99;
  }  
  PrtTim("--- End of arrt ---");
  cout << " -----------  END of ArrayTest (Rc= " << rc << ") --------- " << endl;
  return 0;
}


void tstmtx(int n)
{
  TMatrix<r_4> ma(4,7);
  ma = RegularSequence(10., 2.);
  cout << "\n ::::::  Testing TMatrix - TVector ::::::  " << endl;  
  TMatrix<r_4> sma1 = ma(Range(1,2), Range(2,4));
  if (prtlev > 0) 
    cout << " ----- Matrix ma = \n " << ma << endl;
  sma1 = 3.14;
  if (prtlev > 0) 
    cout << " ----- Matrix ma (Avec sma1=ma(Range(1,2), Range(2,3))=3.14) = \n " << ma << endl;
  TMatrix<r_4> mb(2,3);
  mb = 50.;
  TMatrix<r_4> mc = mb+sma1;
  if (prtlev > 0) 
    cout << " ----- Matrix mc = mb(=50.)+sma1(=3.14) = " << mc << endl;
  TMatrix<r_4> md = (float)100+mc;
  if (prtlev > 0) {
    cout << " ----- Matrix md=100.+mc = " << md << endl;
    cout << " ----- Matrix md-50 = " << md-(float)50. << endl;
  }

  TMatrix<int_4> ima(3,10);
  ima = RegularSequence(0.01,1.);
  if (prtlev > 0) {
    cout << " ----- Matrix ima(3,10) = RegularSequence(0.01,1.) = " << ima << endl;
    cout << " ----- Matrix  ima.Column(3) = \n" << ima.Column(3) << endl;
  }
  TVector<int_4> vcol = ima.Column(7);
  if (prtlev > 0) {
    cout << " ----- Vector  vcol=ima.Column(7) = \n" << vcol << endl;
    cout << " ----- Matrix  ima.Row(0) = \n" << ima.Row(0) << endl;
  }
  TVector<int_4> vlig;
  vlig = ima.Row(1);
  if (prtlev > 0) 
    cout << " ----- Vector  vlig=ima.Row(1) = \n" << vlig << endl;

  TMatrix<int_4> imb(3,10);
  imb = RegularSequence(100.01,10.);
  if (prtlev > 0) 
    cout << " ----- Matrix imb(3,10) = RegularSequence(100.01,10.) = " << imb << endl;
  TMatrix<int_4> imc = ima+imb;
  if (prtlev > 0) 
    cout << " ----- Matrix imc=ima+imb = " << imc << endl;
  TMatrix<int_4> imd(3,10);
  TMatrix<int_4> imd2(3,10);
  imd += 25;
  imd *= 2;
  imd2 = 2;
  imd.MulElt(imd2,imd);
  if (prtlev > 0) 
    cout << " ----- Matrix imd (? = 100) " << imd << endl;
  imb -= imd;
  imc -= imd;
  if (prtlev > 0) 
    cout << " ----- Matrix ((imc-=imd)-(imb-=imd))-ima (== 0 ?) " 
         << (imc-imb)-ima << endl;
  TVector<r_4> va(16, BaseArray::RowVector);
  va = RegularSequence(0.,1.);
  if (prtlev > 0) 
    cout << " ----- Vector va = RegularSequence(0.,1.) = " << va << endl;
  va(Range(2,0,6,2)) *= 10.;
  if (prtlev > 0) 
    cout << " ----- Vector va apres va(Range(2,0,6,2)) *= 10. " << va << endl;
  va(Range(0,5)) += 100.;
  if (prtlev > 0) 
    cout << " ----- Vector va apres va(Range(0,6)) += 100. " << va << endl;

  TMatrix<r_4> esm(3,5);
  EnumeratedSequence ess;
              ess ,
         11 , 12 , 13, 14, 15,  
         21 , 22 , 23, 24, 25,
         31 , 32 , 33, 34, 35  ;  
  esm = ess;
  if (prtlev > 0) 
    cout << "\n TMatrix<r_4> esm(3,5) = EnumeratedSequence(), 11 , 12 , 13, 14, 15 ... " 
         << esm << endl;
  r_4 emin,emax,esum,eprod,esumx2;
  esm.MinMax(emin, emax);  
  esum = esm.Sum();
  eprod = esm.Product();
  esumx2 = esm.SumSq();

  if (prtlev > 0) 
    cout << " Min= " << emin << "  Max= " << emax 
         << " Sum= " << esum << " Prod=" << eprod << " SumSq=" << esumx2 << endl;

  cout << "\n :::::: Testing Matrix operations :::::: " << endl;
  TMatrix<r_4> ra(2,2), ira(2,2);
  /*  ra(0,0) = 2.;   ra(0,1) = 1.;
  ra(1,0) = 4.;   ra(1,1) = 3.;
  */
  ira(0,0) = 1.5;   ira(0,1) = -0.5;
  ira(1,0) = -2.;    ira(1,1) = 1.;
  EnumeratedSequence es;
  es, 2 , 1 , 4 , 3 ;
  ra = es ;
  
  if (prtlev > 0) {
    cout << " ----- Matrix ra = \n " << ra << endl; 
    cout << " ----- Matrix ira = \n " << ira << endl; 
    cout << " ----- Matrix ira*ra = (= Identity ?) " << ira*ra << endl; 
  }

  cout << "\n ::::::: Testing Matrix Inversion ::::::: " << endl;
  SimpleMatrixOperation<r_4> smr;
  TMatrix<r_4> cira = smr.Inverse(ra);
  if (prtlev > 0) {
    cout << " ----- Matrix Inverse(ra) =  " << cira << endl;
    cout << " ----- ira-Inverse(ra) = " << ira-cira << endl;
  }
  Matrix a(n,n);

    
  //  for(int i=0; i<n; i++)
  //    for(int j=0; j<n; j++)  a(j,i) = GauRnd(0., 1.);
  a = RandomSequence(RandomSequence::Gaussian, 0., 2.5);
  if (prtlev > 0) 
    cout << " ----- Matrix a=RandomSequence(RandomSequence::Gaussian, 0., 2.5) = \n " << a << endl;
  Vector x(n);
  x = RegularSequence(1.,3.);
  Vector b = a*x;
  if (prtlev > 0) {
    cout << " ----- Vector x = \n " << x << endl;
    cout << " ----- Vector b = a*x = \n " << b << endl;
  }
  SimpleMatrixOperation<r_8> smo;
  Matrix inva = smo.Inverse(a);
  if (prtlev > 0) {
    cout << " ----- Matrix Inverse(a) = \n " << inva << endl;
    cout << " ----- Matrix a*Inverse(a) = \n " << inva*a << endl;
    cout << " ----- Matrix Inverse(a)*b (=Inv(a)*a*x) = \n " << inva*b << endl;
    cout << " ----- Matrix x-Inverse(a)*b = (=0 ?)\n " << x-inva*b << endl;
  }
}

void tstarr(bool tio)
{
  cout << "\n -----> Testing TArray <---- " << endl;  
  // We create a integer array SizeX=7, SizeY=5
  TArray<int_4> ia(7,5); 
  // We fill it with a sequence of numbers starting at 10., with step = 2.
  ia = RegularSequence(10., 2.);
  cout << " ----- Array IA = \n " << ia << endl;
  // sub array extraction, Range(2,4) : starting position=2 , End=4
  TArray<int_4> ic = ia(Range(2,3),Range(1,3),Range(0));
  cout << " ----- Array IC IA(Range(2,3),Range(1,3)) = \n " << ic << endl;
  // we set the sub-array to zero, this should reflect in the original array
  // sub-arrays share their data with parent array
  ic = 0;
  cout << " ----- Array IC Apres (=0) = \n " << ic << endl;
  cout << " ----- Array IA Apres IC=0 = \n " << ia << endl;

  cout << "  :::: 3 Dim arrays ::::: " << endl;
  TArray<int_4>::SetMaxPrint(1000);
  // Creating 3-dim array (X=8 x Y=7 x Z=2) , filling it with 5
  TArray<int_4> ib(8,7,2);
  ib = 5;
  cout << " ----- Array IB = \n " << ib << endl;
  // Sub array extraction  X from 1 , size 4 - Y from 2 , size 3 , in Z default, from 0, size 1
  // we multiply this sub-array elements by 3
  ib(Range(1,4),Range(2,4), Range(0)) *= 3;
  cout << " -- Array IB , Apres ib(Range(1,4),Range(2,4))*=3 : " << endl;
  cout << ib;
  cout << " -- Array ib(Range(1,4),Range(2,4)) =  : " 
       << ib(Range(1,4),Range(2,4), Range(0)) << endl;

  // Creating a double array X=5 x Y=2
  TArray<r_4> fa(5,2);
  // fill it up with a sequence of 0. to 1.
  fa = RegularSequence(0.,1./(5*2));
  cout << " ------ TArray<r_4> fa(5,2) = \n" << fa << endl;
  // Create a new array from the original array , multiplying it by 2*Pi
  TArray<r_4> fa2 = fa*(float)(2.*M_PI); 
  cout << " ------ TArray<r_4> fa2=fa*2*Pi = \n" << fa2 << endl;
  // Compute sin(fa2) cos(fa2)
  cout << " ------ sin(fa2=fa*2*Pi) = \n" << Sin(fa2) << endl;
  cout << " ------ cos(fa2=fa*2*Pi) = \n" << Cos(fa2) << endl;

  //  Test des fonctions sur complexes 
  TVector<r_4> p_real(10, BaseArray::RowVector), p_imag(10, BaseArray::RowVector);
  p_real = RegularSequence(0., 0.5);
  p_imag = RegularSequence(0., 0.25);
  TVector< complex<r_4> > zvec = ComplexArray(p_real, p_imag);
  cout << " ------ zvec = ComplexArray(p_real, p_imag) -------- "  << endl;
  cout << " :: p_real= " << p_real; 
  cout << " :: p_imag= " << p_imag;
  cout << " :: zvec= " << zvec;
  cout << " :: real(zvec) = " << real(zvec) ;
  cout << " :::: imag(zvec) = " << imag(zvec) ;
  cout << " :::: module2(zvec) = " << module2(zvec) ;
  cout << " :::: module(zvec) = " << module(zvec) ;
  cout << " :::: phase(zvec) = " << phase(zvec) ;
  cout << " ----------------------------------------------------- " << endl;

  if (tio) tstioarr(ia, ib, ic);
}

void tstioarr(TArray<int_4> & ia, TArray<int_4> & ib, TArray<int_4> & ic)
{
  cout << " ------ tstioarr(TArray<int_4> & ia, TArray<int_4> & ib, TArray<int_4> & ic) ---- " << endl;
  {
  cout << " >>>>>> Writing in  arrt.ppf <<<<<<< " << endl;
  POutPersist pos("arrt.ppf");
  // We write  the three arrays in the stream 
  pos << ia << ib << ic;
  cout << " >>>>>> Writing in  arrtn.ppf with names <<<<<<<" << endl;
  POutPersist posn("arrtn.ppf");
  string tag = "ArrIA";
  posn.PutObject(ia, tag);
  tag = "ArrIB";
  posn.PutObject(ib, tag);
  tag = "ArrIC";
  posn.PutObject(ic, tag);
  }

  {
  cout << " >>>>>>> Reading from  arrt.ppf  <<<<< " << endl;
  PInPersist pis("arrt.ppf");
  TArray<int_4> iaa, ibb, icc; 
  // We read the three arrays from the stream 
  pis >> iaa >> ibb >> icc;
  cout << " ----- Array IAA = \n " << iaa << endl;
  cout << " ----- Array IBB = \n " << ibb << endl;
  cout << " ----- Array ICC = \n " << icc << endl;
  icc = 12;
  cout << " ----- Array ICC (=12) = \n " << icc << endl;
  cout << " ----- Array IAA (ICC=12) = \n " << iaa << endl;
  }

  {
  cout << " >>>>>>> Reading from  arrtn.ppf  <<<<< " << endl;
  PInPersist pis("arrtn.ppf");
  TArray<int_4> iaa, ibb, icc; 
  // We read the three arrays from the stream 
  string tag = "ArrIC";
  pis.GetObject(icc, tag);
  tag = "ArrIB";
  pis.GetObject(ibb, tag);
  tag = "ArrIA";
  pis.GetObject(iaa, tag);

  cout << " ----- Array IAAA = \n " << iaa << endl;
  cout << " ----- Array IBBB = \n " << ibb << endl;
  cout << " ----- Array ICCC = \n " << icc << endl;
  icc = 68;
  cout << " ----- Array ICCC (=12) = \n " << icc << endl;
  cout << " ----- Array IAAA (ICC=12) = \n " << iaa << endl;
  }

}


void tstasciiio()
{

cout << "\n -----> Testing EnumeratedSequence/TArray ASCII I/O <---- " << endl; 

EnumeratedSequence ess;
int nbad;
ess.Append("56.5 (1.,-1.) 4 8 16", nbad);
cout << " --- ess.Append(\"56.5 (1.,-1.) 4 8 16\", nbad) :  ess=" << endl;
cout << ess;

EnumeratedSequence ess2;
ess2.Append("56.5 (-3.14, 12) 'Hello' 4. 8.565432 12 32 3.1415e-6", nbad);
cout << " --- ess2.Append(...) ess2=" << endl; 
cout << ess2;
cout << " --- ess2.Append(...) (2) ess2=" << endl; 
ess2.Append("(2.7, -8) Bonjour ByeBye 'Comment ca va ?' 64 64.7 128 1.28 25.6e3", 
	    nbad);
cout << ess2;

Matrix A(4,6);
A = RegularSequence(0.5, 0.2);
cout << " Matrix A(4,6); A=RegularSequence(0.5, 0.2) : " << endl;
cout << A ;

char * flnm = "arr_ascii.txt";
{
cout << " Writing A to file " << flnm << endl;
ofstream ofs(flnm);
A.WriteASCII(ofs);
}

Matrix B;
sa_size_t nr, nc;
{
ifstream ifs(flnm);
B.ReadASCII(ifs,nr,nc);
}

cout << " Matrix B; B.ReadASCII() from file " << flnm << endl;
cout << B ;

Vector C;
{
ifstream ifs(flnm);
ifs >> C;
}
cout << " Vector C; ifstream( " << flnm << ") >> C " << endl;
cout << " C.Transpose()= " << C.Transpose() ;

TVector<int_4> IC;
{
ifstream ifs(flnm);
IC.ReadASCII(ifs,nr,nc);
}
cout << " TVector<int_4> IC; IC.ReadASCII() from file " << flnm << endl;
cout << " IC.Transpose()= " << IC.Transpose() ;

TArray<r_4> D(4,6,2);
{
ifstream ifs(flnm);
D.ReadASCII(ifs,nr,nc);
}
cout << " TArray<r_4> D(4,6,2); D.ReadASCII() from file " << flnm << endl;
cout << D ;

}

// ----- Fonction de verification des fonctions de base de TArray 
// Retourne 0 si OK 
int checktarr()
{
  int rc = 0;

  int NX = 17;
  int NY = 8;
  int NZ = 6;
  int NT = 5;
  if (prtlev > 0)  
    cout << " CheckArry: Starting Test-1 EltAcc/SubArray <int> " << endl;
  TArray<int> IA(NX, NY, NZ, NT);
  for(int ix=0; ix<NX; ix++)
    for(int iy=0; iy<NY; iy++)
      for(int iz=0; iz<NZ; iz++)
	for(int it=0; it<NT; it++)
	  IA(ix, iy, iz, it) = it*10000+iz*1000+iy*100+ix;

  TArray<int> IB(NX, NY);
  TMatrix<int> IMXB(IB);
  for(int ir=0; ir<IMXB.NRows(); ir++) {
    IMXB.Row(ir) = RegularSequence(0.1);  
    IMXB.Row(ir) += ir*100;
  }

  int nerr = 0;
  int min=99, max=99;
  for(int iz=0; iz<NZ; iz++)
    for(int it=0; it<NT; it++)  {
      TArray<int> IC = IB+(it*10000+iz*1000);
      TArray<int> IAI = IA(Range::all(), Range::all(), Range(iz, iz), Range(it, it));
      // Fait par SubArray      IAI.CompactTrailingDim();
      if (prtlev > 2) {
	cout << " CheckArry/Debug: Test-1 iz=" << iz << " it=" << it << endl;
	cout << "IAI= \n " << IAI << " IC= \n " << IC << endl;
      }
      TArray<int> ID = IAI - IC;
      ID.MinMax(min, max);
      if ((min!=0) || (max!=0)) nerr++;
    }
  
  if (nerr == 0)
    cout << "(1) CheckTArray <int>  EltAcc, SubArray OK " << endl; 
  else {
    cout << "(1) CheckTArray <int>  EltAcc, SubArray ERROR -> Rc += 1 !!! NErr=" 
	 << nerr << endl; 
    rc += 1;
  }

  if (prtlev > 0)  
    cout << " CheckArry: Starting Test-2 SubArray/Range <int> " << endl;  nerr = 0;
  {
  TArray<int> ID, IAI;
  IAI = IA(Range::all(), Range::all(), Range::last(), Range::first());
  // Fait par SubArray  IAI.CompactTrailingDim();
  // cout << "---DBG- IAI= \n " << IAI << " IB= \n " << IB << endl;
  ID = IAI - IB;
  ID -= ((NZ-1)*1000);
  // cout << "---DBG- ID= \n" << ID ;
  ID.MinMax(min, max);
  if ((min!=0) || (max!=0)) nerr++;
  }
  {
  TArray<int> ID, IAI;
  IAI = IA(Range::all(), Range::all(), Range::last(), Range::last());
  // Fait par SubArray  IAI.CompactTrailingDim();
  //  cout << "---DBG- IAI= \n " << IAI << " IB= \n " << IB << endl;
  ID = IAI - IB;
  ID -= ((NZ-1)*1000+(NT-1)*10000);
  // cout << "---DBG- ID= \n" << ID ;
  ID.MinMax(min, max);
  if ((min!=0) || (max!=0)) nerr++;
  }
  if (nerr == 0)
    cout << "(2) CheckTArray <int>  SubArray/Range::last() OK " << endl; 
  else {
    cout << "(2) CheckTArray <int>  SubArray/Range::last() ERROR  Rc += 2 !!! NErr=" 
	 << nerr << endl; 
    rc += 2;
  }

  if (prtlev > 0)  
    cout << " CheckArry: Starting Test-3 SubArray(SubArray) <int> " << endl;  nerr = 0;
  TArray<int> IAE = IA(Range::all(), Range::all(), 
		       //	      		       Range(1, 0, 3, 2), Range::last() );
  		       Range(Range::lastIndex(), Range::lastIndex(), 3, 2), Range::last() );
  // Fait par SubArray  IAE.CompactTrailingDim();

  TArray<int> IBE(NX/3, NY/2);  
  for(int ix=0; ix<NX/3; ix++)
    for(int iy=0; iy<NY/2; iy++)
      IBE(ix, iy, 0) = 3*ix+(2*100*iy);

  nerr = 0;
  for(int iz=0; iz<IAE.SizeZ(); iz++) {
    TArray<int> IAI = IAE(Range::all(), Range::all(), Range(iz, iz));
    // Fait par SubArray  IAI.CompactTrailingDim();
    //    cout << "---DBG- iz= " << iz << " IAI= \n " << IAI << " IB= \n " << IB << endl;
    TArray<int> ID = IAI-IB;
    ID -= ((iz*2+1)*1000+(NT-1)*10000);
    ID.MinMax(min, max);
    if ((min!=0) || (max!=0)) nerr++;
  }

  if (nerr == 0)
    cout << "(3) CheckTArray <int>  SubArray(SubArray) OK " << endl; 
  else {
    cout << "(3) CheckTArray <int>  SubArray(SubArray)  ERROR  Rc += 4 !!! NErr=" 
	 << nerr << endl; 
    rc += 4;
  }

  return rc;
}