source: Sophya/trunk/Poubelle/DPC:FitsIOServer/NTools/perandom.cc@ 927

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

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

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

//++
FunRan::FunRan(double *tab, int nBin)
//
//      Createur.
//--
: Histo(0, (float)(nBin), nBin)
{
  (*this)(0) = tab[0];
  for(int 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 j=0; j<nBin; j++)
    (*this)(j) /= (*this)(nBin-1);
  END_CONSTRUCTOR
}

FunRan::FunRan(double *tab, int nBin, float xMin, float xMax)
: Histo(xMin, xMax, nBin)
{
  (*this)(0) = tab[0];
  for(int 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 j=0; j<nBin; j++)
    (*this)(j) /= (*this)(nBin-1);
  END_CONSTRUCTOR
}

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

  return iBin;
}

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

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

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

                   

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

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

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

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

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