source: Sophya/trunk/SophyaLib/NTools/perandom.cc@ 2615

#include "sopnamsp.h"
#include "machdefs.h"
#include "pexceptions.h"
#include "perandom.h"
#include "pemath.h"
#include <iostream>

////////////////////////////////////////////////////////////////////////////
//++
// Class        FunRan
// Lib  Outils++ 
// include      perandom.h
//
//      Tirage aleatoire sur un histogramme 1D.
//--

//++
FunRan::FunRan(FunRan::Func f, r_8 xMin, r_8 xMax, int_4 nBin)
//
//      Createur.
//--
: Histo(xMin, xMax, nBin)
{
  (*this)(0) = f(BinLowEdge(0));
  for(int_4 i=1; i<nBin; i++)
    (*this)(i) = (*this)(i-1) + f(BinLowEdge(i));
    
  for(int_4 j=0; j<nBin; j++)
    (*this)(j) /= (*this)(nBin-1);
}

//++
FunRan::FunRan(r_8 *tab, int_4 nBin)
//
//      Createur.
//--
: Histo(0, (r_8)(nBin), nBin)
{
  (*this)(0) = tab[0];
  for(int_4 i=1; i<nBin; i++)
    (*this)(i) = (*this)(i-1) + tab[i];

  if ((*this)(nBin-1) == 0) {
    cerr << "FunRan::FunRan : sum(prob) = 0" << endl;
    exit(-1);
  }

  for(int_4 j=0; j<nBin; j++)
    (*this)(j) /= (*this)(nBin-1);
}

FunRan::FunRan(r_8 *tab, int_4 nBin, r_8 xMin, r_8 xMax)
: Histo(xMin, xMax, nBin)
{
  (*this)(0) = tab[0];
  for(int_4 i=1; i<nBin; i++)
    (*this)(i) = (*this)(i-1) + tab[i];

  if ((*this)(nBin-1) == 0) {
    cerr << "FunRan::FunRan : sum(prob) = 0" << endl;
    exit(-1);
  }

  for(int_4 j=0; j<nBin; j++)
    (*this)(j) /= (*this)(nBin-1);
}

//++
int_4 FunRan::BinRandom()
//
//      Tirage avec retour du numero de bin.
//--
{
  r_8 z=drand01();
  if (z <= 0) return 0;
  if (z >= 1) return mBins-1;
  
  // recherche du premier bin plus grand que z
  int_4 iBin = 0;
  for (; iBin<mBins; iBin++)
    if (z < (*this)(iBin)) break;

  return iBin;
}

//++
r_8 FunRan::Random()
//
//      Tirage avec retour abscisse du bin interpole.
//--
{
  r_8 z=drand01();
  if (z <= 0) return mMin;
  if (z >= 1) return mMax;
  // cas z <= tab[0] 
  if (z <= (*this)(0)) {
    r_8 t = mMin + binWidth/(*this)(0) * z;
    return t;
  }

  // recherche du premier bin plus grand que z
  int_4 iBin = 0;
  for (; iBin<mBins; iBin++)
    if (z < (*this)(iBin)) break;

  // interpolation pour trouver la valeur du tirage aleatoire
  r_8 t1 = (*this)(iBin-1);
  r_8 x1 = BinLowEdge(iBin-1);
  r_8 t2 = (*this)(iBin);
  r_8 x2 = x1 + binWidth;
  r_8 t = x1 + (x2-x1) / (t2-t1) * (z-t1);
  if (t < mMin) t = mMin;
  if (t > mMax) t = mMax;
  return(t);
}

                   

////////////////////////////////////////////////////////////////////////////
//++
// Class        FunRan2D
// Lib  Outils++ 
// include      perandom.h
//
//      Tirage aleatoire sur un histogramme 2D.
//--

//++
FunRan2D::FunRan2D(r_8 *tab, int_4 nBinX, int_4 nBinY)
//
//      Createur.
//--
{
  // Tirage en X, somme sur les Y.
   r_8* tabX = new r_8[nBinX];
   for (int_4 i=0; i<nBinX; i++) {
     tabX[i] = 0;
     for (int_4 j=0; j<nBinY; j++) {
       tabX[i] += tab[i*nBinY +j];
     }
   }
   ranX = new FunRan(tabX, nBinX);
   delete[] tabX;
   
   ranY = new(FunRan*[nBinX]);
   
   for (int_4 k=0; k<nBinX; k++)
      ranY[k] = new FunRan(tab + nBinY*k, nBinY);
   
   nx = nBinX;
}

//++
FunRan2D::FunRan2D(r_8 **tab, int_4 nBinX, int_4 nBinY)
//
//      Createur.
//--
{
  // Tirage en X, somme sur les Y.
   r_8* tabX = new r_8[nBinX];
   for (int_4 i=0; i<nBinX; i++) {
     tabX[i] = 0;
     for (int_4 j=0; j<nBinY; j++) {
       tabX[i] += tab[i][j];
     }
   }
   ranX = new FunRan(tabX, nBinX);
   
   ranY = new(FunRan*[nBinX]);
   
   for (int_4 k=0; k<nBinX; k++)
    if (tabX[k] != 0) 
      ranY[k] = new FunRan(tab[k], nBinY);
    else ranY[k] = NULL;
   
   delete[] tabX;
   nx = nBinX;
}

FunRan2D::~FunRan2D()
{
  for (int_4 i=nx-1; i>=0; i--)
    delete ranY[i];
    
  delete[] ranY;
  
  delete ranX;
}

//++
void FunRan2D::BinRandom(int_4& x, int_4& y)
//
//      Tirage avec retour du numeros de bin.
//--
{
  x = ranX->BinRandom();
  //  FAILNIL(ranY[x]);  Ne compile pas $CHECK$ Reza 22/04/99
  y = ranY[x]->BinRandom();
}

//++
void FunRan2D::Random(r_8& x, r_8& y)
//
//      Tirage avec retour abscisse et ordonnee
//      du bin interpole.
//--
{
  x = ranX->Random();
  int_4 i = int_4(ceil(x));
  //  FAILNIL(ranY[i]);  Ne compile pas $CHECK$ Reza 22/04/99
  y = ranY[i]->Random();
}