source: Sophya/trunk/SophyaLib/BaseTools/stsrand.cc@ 3401

#include "stsrand.h"
#include "thsafeop.h"
#include "srandgen.h"
#include "fiondblock.h"
#include <math.h> 

namespace SOPHYA {

/*!
   \class RandomGenerator
   \ingroup BaseTools
   This class is a thread-safe random number generator.
   Its PPF handler can be used to save the complete state of the class and the underlying
   random number generator used.

   \sa SOPHYA::ObjFileIO<RandomGenerator>

*/

// Objet statique global pour gestion de lock entre threads
static ThSafeOp* ths_rand = NULL;

RandomGenerator::RandomGenerator(size_t n, bool tsafe) 
{
  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
  if (tsafe) {   // thread-safe
    fg_nothrsafe = false;
    if (n < 1) n = 1024;
    rseq_.ReSize(n, false);
    idx_ = n;
  }
  else {   // NOT thread-safe
    fg_nothrsafe = true;
    idx_ = 1;
  }
}

RandomGenerator::RandomGenerator(RandomGenerator const & rg)
{
  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
  if (!rg.fg_nothrsafe) { // thread-safe
    fg_nothrsafe = false;
    rseq_.ReSize(rg.rseq_.Size(), false);
    idx_ = rseq_.Size();
  }
  else {   // NOT thread-safe
    fg_nothrsafe = true;
    idx_ = 1;
  }
}
 

RandomGenerator::~RandomGenerator(void) 
{
  // rien a faire 
}

void RandomGenerator::AutoInit(int lp)
{
  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
  ths_rand->lock();
  if(lp) cout << "RandomGenerator::AutoInit() : Calling Auto_Ini_Ranf() ..." << endl;
  Auto_Ini_Ranf(lp);  // Faut-il faire copier/coller du code Auto_Ini_Ranf() ici ? 
  ths_rand->unlock();
}

void RandomGenerator::Init(long seed_val, int lp)
{
  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
  if(lp) cout << "RandomGenerator::Init(long seed=" << seed_val << ")" << endl;
  ths_rand->lock();
  srand48(seed_val);
  ths_rand->unlock();
  return;
}

void RandomGenerator::Init(unsigned short seed_16v[3], int lp)
{
  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
  if(lp) cout << "RandomGenerator::Init(u_short seed_16v[3]=" << seed_16v[0]
              << "," << seed_16v[1] << "," << seed_16v[2] << ")" << endl;
  ths_rand->lock();
  Init_P(seed_16v);
  ths_rand->unlock();
}

void RandomGenerator::GetSeed(unsigned short seed_16v[3], int lp)
{
  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
  ths_rand->lock();
  GetSeed_P(seed_16v);
  ths_rand->unlock();
  if(lp) cout << "RandomGenerator::GetSeed(u_short seed_16v[3]=" << seed_16v[0]
              << "," << seed_16v[1] << "," << seed_16v[2] << ")" << endl;
  return;
}

void RandomGenerator::Init_P(unsigned short seed_16v[3])
{
  seed48(seed_16v);
}

void RandomGenerator::GetSeed_P(unsigned short seed_16v[3])
{
  unsigned short seed[3] = {0,0,0};
  unsigned short *p;
  p = seed48(seed);
  memcpy(seed_16v,p,3*sizeof(unsigned short));
  /* on re-initialise a ce qui etait avant */
  seed48(seed_16v);
  return;
}

r_8 RandomGenerator::Gaussian() 
{
  r_8 A=Next(); 
  while (A==0.) A=Next(); 
  return sqrt(-2.*log(A))*cos(2.*M_PI*Next());
}


uint_8 RandomGenerator::Poisson(double mu,double mumax)
{
  double pp,ppi,x;

  if((mumax>0.)&&(mu>=mumax)) {
    pp = sqrt(mu);
    while( (x=pp*Gaussian()) < -mu );
    return (uint_8)(mu+x+0.5);
  }
  else {
    uint_8 n;
    ppi = pp = exp(-mu);
    x = Next();
    n = 0;
    while (x > ppi) {
      n++;
      pp = mu*pp/(double)n;
      ppi += pp;
    }
    return n;
  }
  return 0;  // pas necessaire ?
}

void RandomGenerator::GenSeq(void)
{
  ths_rand->lock();
  for(size_t k=0; k<rseq_.Size(); k++)  rseq_(k) = drand48();
  ths_rand->unlock();
  idx_ = 0;
}

//----------------------------------------------------------
// Classe pour la gestion de persistance
// ObjFileIO<RandomGenerator>
//----------------------------------------------------------

/* --Methode-- */
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void        ObjFileIO<RandomGenerator>::WriteSelf(POutPersist& s) const
{
  if (dobj == NULL)
    throw NullPtrError("ObjFileIO<RandomGenerator>::WriteSelf() dobj=NULL");
  ths_rand->lock();  // thread-safety
  uint_4 itab[6];
  //itab : [0]: version, [1,2,3] = srand48 state/seed  , [4,5] = reserved for future use
  itab[0] = 1;
  // On recupere et on ecrit ds le PPF l'etat du generateur aleatoire
  unsigned short seed_16v[3];
  dobj->GetSeed_P(seed_16v);
  for(int i=0; i<3; i++)  itab[i+1] = seed_16v[i];
  itab[4] = 0;
  s.Put(itab, 6);
  uint_8 sz = dobj->rseq_.Size();
  s.Put(sz);  // Taille du tableau intermediaire
  uint_8 ix = dobj->idx_;
  s.Put(ix);  // valeur de l'index 

  if (dobj->rseq_.Size() > 0)  s << dobj->rseq_;  // On ecrit le tableau (NDataBlock) si necessaire
  ths_rand->unlock(); // thread-safety
  return;
}

/* --Methode-- */
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void        ObjFileIO<RandomGenerator>::ReadSelf(PInPersist& s)
{
  uint_4 itab[6];
  //itab : [0]: version, [1,2,3] = srand48 state/seed  , [4] = reserved for future use
  s.Get(itab, 6);
  uint_8 sz,ix;
  s.Get(sz);  // Taille du tableau intermediaire
  s.Get(ix);  // Taille du tableau intermediaire

  if (dobj == NULL) dobj = new RandomGenerator(sz, (sz>0)?true:false);
  dobj->idx_ = ix;
  if (sz > 0) { 
    s >> dobj->rseq_;  // On lit le tableau (NDataBlock) si necessaire
    dobj->fg_nothrsafe = false;
  }
  else {  // Objet lu est NON thread-safe, taille_tableau rseq_ = 0
    dobj->fg_nothrsafe = true;
    if (dobj->rseq_.Size() > 0)  dobj->rseq_.Dealloc(); 
  }
  // On initialise l'etat du generateur aleatoire avec les valeurs lues
  unsigned short seed_16v[3];
  dobj->GetSeed_P(seed_16v);
  for(int i=0; i<3; i++)  seed_16v[i] = itab[i+1];
  dobj->Init(seed_16v, 0);
  return;
}

// ---------------------------------------------------------
#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template ObjFileIO<RandomGenerator>
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class ObjFileIO<RandomGenerator>;
#endif
// ---------------------------------------------------------

}  /* namespace SOPHYA */