// Test des operations de tirages aleatoires thread safe/nosafe
#include "sopnamsp.h"
#include "zthread.h"

#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>
#include <string>
#include <iostream>

#include "tvector.h"
#include "matharr.h"
#include "fioarr.h"
#include "tarrinit.h"
#include "srandgen.h"
#include "nbtri.h"
#include "timing.h"

#include "randr48.h"

long Check_Redondances(TVector<r_8>& v);

//---------------------------------------------------------------------
// Classe heritant de ZThread, pour remplir un vecteur avec des aleatoires
class MTRndTest : public ZThread {
public:
  MTRndTest(TVector<r_8>& v,bool thsafe=true,int nbuff=1024);
  virtual void run(); 
protected:
  string nom_; int inum_;
  TVector<r_8> vv_;
  bool thsafe_;
  int nbuff_;
};

static int mtrandId = 0;  // Pour donner un identificateur a chaque thread

MTRndTest::MTRndTest(Vector& v,bool thsafe,int nbuff)
  : vv_(v,true) , thsafe_(thsafe), nbuff_(nbuff)
{
  if(nbuff_<=0) nbuff_ = 1024;
  char buff[32];
  inum_ = mtrandId;
  sprintf(buff,"MTRndTest-Id=%d",mtrandId);
  mtrandId++;
  nom_ = buff;
  cout<<"Thread MTRndTest("<<nom_<<" ) Created ... with threadsafe="<<thsafe_
      <<", nbuff="<<nbuff_<<endl;
}


// Le travail fait dans chaque thread (appelle par start())
void MTRndTest::run()
{
  cout<<"MTRndTest::run() - Nom= "<<nom_<<" vv.Size()= "<<vv_.Size()
      <<" with threadsafe="<<thsafe_<<endl;
  if(thsafe_) {  // On appele ThSDR48RandGen Thread-Safe
    ThSDR48RandGen rgen(nbuff_,true);
    for(sa_size_t k=0; k<vv_.Size(); k++) vv_(k) = rgen.Flat01();
  } else {  // On appele ThSDR48RandGen NO Thread-Safe
    // ATTENTION: le comportement est different selon qu'un thread
    // se termine avant ou apres que le suivant commence:
    // Pour de petits VSZ le run() du thread "n" se termine
    // avant que le run() du thread "n+1" ne commence:
    // cela donne un comportement "tread_safe" APPARENT mais
    // dans la realite drand48() N'EST PAS thread_safe !!!!
    ThSDR48RandGen rgen(0,false);
    for(sa_size_t k=0; k<vv_.Size(); k++) vv_(k) = rgen.Flat01();
  }
  cout<<"....MTRndTest::run() - "<<nom_<<" ENDED !!!!!!!!!"<<endl;
}

//---------------------------------------------------------------------
int main(int narg, char *arg[])

{
  SophyaInit();

  if(narg<3) {cout<<"usage: tmtrnd3 ntir nthreads [nbuff]"<<endl; return -1;}
  sa_size_t VSZ = atoi(arg[1]);
  int NTH = atoi(arg[2]);
  int NBUFF = (narg>3) ? atoi(arg[3]): 64;
  if(NTH<=0) NTH = 1;
  cout<<"tmtrnd/arguments, VSZ= "<<VSZ<<" NTH= "<<NTH<<" NBUFF="<<NBUFF<<endl;
  // Il ne faut pas que les buffers de ThSDR48RandGen soient remplis
  // au dela de ce qui est necessaire si on veut faire
  // la comparaison avec le mono thread
  int lmod = NTH*NBUFF;
  if(VSZ<lmod) VSZ = lmod;
  if(VSZ%lmod!=0) VSZ = int(VSZ/(double)lmod+0.5)*lmod;
  cout<<"tmtrnd/starting, lmod="<<lmod<<" NTH= "<<NTH<<" -> changed VSZ= "<<VSZ<<endl;

  int rc = 0;
  sa_size_t csz = VSZ / NTH;
  sa_size_t first, last;

  try {

    //---
    //--- Init et set de l'aleatoire
    //---
    ThSDR48RandGen rgen_nothsf(0,false);
    unsigned short seed_16v[3];
    {
    cout<<"\n>>>>>>> Init de l'aleatoire"<<endl;
    rgen_nothsf.AutoInit(5);
    rgen_nothsf.GetSeed(seed_16v);
    }

    InitTim();

    //---
    //--- remplissage en mono-thread directement avec drand48
    //---
    TVector<r_8> DATA0(VSZ);
    {
    cout<<"\n>>>>>>> Remplissage en mono-thread avec drand48()"<<endl;
    rgen_nothsf.SetSeed(seed_16v);
    for(int i=0;i<VSZ;i++) DATA0(i) = drand48();
    PrtTim(" ");
    }

    //---
    //--- Remplissage par threads avec ThSDR48RandGen thread safe
    //---
    TVector<r_8> DATA1(VSZ);
    {
    cout<<"\n>>>>>>> Remplissage par threads avec ThSDR48RandGen thread safe"<<endl;
    rgen_nothsf.SetSeed(seed_16v);
    vector<MTRndTest *> vth1;
    cout<<"...tmtrnd/creating threads "<<endl;
    for(int kt=0; kt<NTH; kt++) {
      first = kt*csz;
      last = (kt == NTH-1) ? VSZ-1 : first+csz-1;
      TVector<r_8> sv = DATA1(Range(first, last));
      vth1.push_back(new MTRndTest(sv,true,NBUFF) );
    }
    cout<<"...tmtrnd/starting threads"<<endl;
    for(int kt=0; kt<NTH; kt++)  vth1[kt]->start();
    cout<<"...tmtrnd/waiting for all  threads to finish"<<endl;
    sleep(2);
    for(int kt=0; kt<NTH; kt++)  vth1[kt]->join();
    cout<<"...tmtrnd/deleting thread objects"<<endl;
    for(int kt=0; kt<NTH; kt++)  delete vth1[kt];
    PrtTim(" ");
    }

    //---
    //--- Remplissage par threads avec ThSDR48RandGen SANS thread safe
    //---
    TVector<r_8> DATA2(VSZ);
    {
    cout<<"\n>>>>>>> Remplissage par threads avec ThSDR48RandGen SANS thread safe"<<endl;
    rgen_nothsf.SetSeed(seed_16v);
    vector<MTRndTest *> vth2;
    cout << "......tmtrnd/creating threads " << endl;
    for(int kt=0; kt<NTH; kt++) {
      first = kt*csz;
      last = (kt == NTH-1) ? VSZ-1 : first+csz-1;
      TVector<r_8> sv = DATA2(Range(first, last));
      vth2.push_back(new MTRndTest(sv,false) );
    }
    cout<<"......tmtrnd/starting threads"<<endl;
    for(int kt=0; kt<NTH; kt++)  vth2[kt]->start();
    cout<<"......tmtrnd/waiting for all  threads to finish"<<endl;
    sleep(2);
    for(int kt=0; kt<NTH; kt++)  vth2[kt]->join();
    cout<<"......tmtrnd/deleting thread objects "<<endl;
    for(int kt=0; kt<NTH; kt++)  delete vth2[kt];
    PrtTim(" ");
    }

    //---
    //--- remplissage en mono-thread avec ThSDR48RandGen
    //---
    TVector<r_8> DATA0R(VSZ);
    {
    cout<<"\n>>>>>>> Remplissage en mono-thread avec ThSDR48RandGen"<<endl;
    rgen_nothsf.SetSeed(seed_16v);
    for(int i=0;i<VSZ;i++) DATA0R(i) = rgen_nothsf.Flat01();
    PrtTim(" ");
    }

    //---
    //--- checking
    //---
    cout<<"\n>>>>>>> Sorting"<<endl;
    qsort(DATA0.Data(),(size_t)VSZ,sizeof(double),qSort_Dble);
    qsort(DATA1.Data(),(size_t)VSZ,sizeof(double),qSort_Dble);
    qsort(DATA2.Data(),(size_t)VSZ,sizeof(double),qSort_Dble);
    qsort(DATA0R.Data(),(size_t)VSZ,sizeof(double),qSort_Dble);
    PrtTim(" ");

    cout<<"\n>>>>>>> Checking differences"<<endl;
    double vmax00=0., vmax01=0., vmax02=0., vmax12=0.;
    for(int i=0;i<VSZ;i++) {
      double v;
      v = fabs(DATA0R(i)-DATA0(i));
      if(v>vmax00) vmax00 = v;
      v = fabs(DATA1(i)-DATA0(i));
      if(v>vmax01) vmax01 = v;
      v = fabs(DATA2(i)-DATA0(i));
      if(v>vmax02) vmax02 = v;
      v = fabs(DATA2(i)-DATA1(i));
      if(v>vmax12) vmax12 = v;
    }
    cout<<"vmax00 = "<<vmax00<<endl
        <<"vmax01 = "<<vmax01<<endl
        <<"vmax02 = "<<vmax02<<endl
        <<"vmax12 = "<<vmax12<<endl;

    
    cout<<"\n>>>>>>> Checking redondances"<<endl;
    cout<<"DATA0:"<<endl; Check_Redondances(DATA0);
    cout<<"DATA1:"<<endl; Check_Redondances(DATA1);
    cout<<"DATA2:"<<endl; Check_Redondances(DATA2);
    cout<<"DATA0R:"<<endl; Check_Redondances(DATA0R);

    //--- persisting
    cout<<"\n>>>>>>> Writing to ppf"<<endl;
    string tag;
    POutPersist pos("tmtrnd3.ppf");
    tag = "d0r"; pos.PutObject(DATA0R,tag);
    tag = "d0"; pos.PutObject(DATA0,tag);
    tag = "d1"; pos.PutObject(DATA1,tag);
    tag = "d2"; pos.PutObject(DATA2,tag);

    double eps = 1.e-25;
    if(vmax00>eps || vmax01>eps) rc = 1;

  }
  catch (PThrowable exc) {
    cerr<<"zthr: catched Exception "<<exc.Msg()<<endl;
    rc = 77;
  }  
  catch (...) {
    cerr<<"catched unknown (...) exception (lpk.cc)"<<endl; 
    rc = 78; 
  } 

  cout<<"\n>>>>>>>  tmtrnd3/END"<<endl;

  return(rc);
}

//---------------------------------------------------------------------
long Check_Redondances(TVector<r_8>& v)
// Recherche des aleatoire tires 2 fois
// Le vecteur "v" doit deja etre trie par valeurs CROISSANTES
{
 int k=0;
 double eps=1.e-25;
 TVector<r_8> dd(v,false);
 for(int i=dd.Size()-1;i>=1;i--) dd(i) -= dd(i-1);
 dd(0) = 0.;
 qsort(dd.Data(),(size_t)dd.Size(),sizeof(double),qSort_Dble);
 for(int i=1;i<dd.Size();i++) if(fabs(dd(i))<eps) k++;
 cout<<"Number of redondances pour v<"<<eps<<" : k="<<k<<endl;
 return k;
}

/*
delobjs *
openppf tmtrnd3.ppf

#------ les tableaux sont ranges par valeurs croissantes
disp d0
disp d0r "same"
disp d1 "same red"
disp d2 "same blue"


#------
newh1d hd0 0. 1. 1000
n/proj hd0 d0.val
newh1d hd1 0. 1. 1000
n/proj hd1 d1.val
newh1d hd2 0. 1. 1000
n/proj hd2 d2.val

disp hd0
disp hd1 "same red"
disp hd2 "same blue"

#------ tableaux des differences / a la reference
del d00 d10 d20 d21
c++exec \
TVector<r_8> d00(d0r,false); d00 -= d0; \
TVector<r_8> d10(d1,false); d10 -= d0; \
TVector<r_8> d20(d2,false); d20 -= d0; \
TVector<r_8> d21(d2,false); d21 -= d1; \
KeepObj(d00); KeepObj(d10); KeepObj(d20); KeepObj(d21); \
cout<<"-- Fin"<<endl;

# Doivent avoir min=max=0
echo d00 min =  ${d00.min} , max =  ${d00.max}
echo d10 min =  ${d10.min} , max =  ${d10.max}
# Doivent avoir min!=0 et max!=0
echo d20 min =  ${d20.min} , max =  ${d20.max}
echo d21 min =  ${d21.min} , max =  ${d21.max}

disp d00
disp d10

disp d20
disp d21 "same red"

#------ tableaux des redondances
del dd1 dd2
c++exec \
TVector<r_8> dd1(d1,false); \
for(int i=dd1.Size()-1;i>=1;i--) dd1(i) -= dd1(i-1); dd1(0) = 0.; \
TVector<r_8> dd2(d2,false); \
for(int i=dd2.Size()-1;i>=1;i--) dd2(i) -= dd2(i-1); dd2(0) = 0.; \
qsort(dd1.Data(),(size_t)dd1.Size(),sizeof(double),qSort_Dble); \
qsort(dd2.Data(),(size_t)dd2.Size(),sizeof(double),qSort_Dble); \
double eps=1.e-25; \
int k1=0; for(int i=1;i<dd1.Size();i++) if(fabs(dd1(i))<eps) k1++; \
int k2=0; for(int i=1;i<dd2.Size();i++) if(fabs(dd2(i))<eps) k2++; \
cout<<"Number of v<"<<eps<<" : k1="<<k1<<" , k2="<<k2<<endl; \
KeepObj(dd1); KeepObj(dd2); \
cout<<"-- Fin"<<endl;

# Regarder pour les petites separations
# dd1 ne doit pas avoir de redondance
disp dd1
# dd2 doit avoir des redondances
disp dd2

 */