// @(#)root/unuran:$Id$
// Authors: L. Moneta,  Dec 2010

/**********************************************************************
 *                                                                    *
 * Copyright (c) 2010  LCG ROOT Math Team, CERN/PH-SFT                *
 *                                                                    *
 *                                                                    *
 **********************************************************************/

// Implementation file for class TFoamSampler

#include "TFoamSampler.h"
#include "Math/DistSamplerOptions.h"

#include "TFoam.h"
#include "TFoamIntegrand.h"
#include "Math/OneDimFunctionAdapter.h"
#include "Math/IOptions.h"
#include "Fit/DataRange.h"

#include "TRandom.h"
#include "TError.h"

#include "TF1.h"
#include <cassert>
#include <cmath>

class FoamDistribution : public TFoamIntegrand {

public:

   FoamDistribution(const ROOT::Math::IMultiGenFunction & f, const ROOT::Fit::DataRange & range) :
      fFunc(f),
      fX(std::vector<double>(f.NDim() ) ),
      fMinX(std::vector<double>(f.NDim() ) ),
      fDeltaX(std::vector<double>(f.NDim() ) )
   {
      assert(f.NDim() == range.NDim() );
      std::vector<double> xmax(f.NDim() );
      for (unsigned int i = 0; i < range.NDim(); ++i) {
         if (range.Size(i) == 0)
            Error("FoamDistribution","Range is not set for coordinate dim %d",i);
         else if (range.Size(i)>1)
            Warning("FoamDistribution","Using only first range in coordinate dim %d",i);

         std::pair<double,double> r = range(i);
         fMinX[i] = r.first;
         fDeltaX[i] = r.second - r.first;
      }
   }
   // in principle function does not need to be cloned

   virtual double Density(int ndim, double * x) {
      assert(ndim == (int) fFunc.NDim() );
      for (int i = 0; i < ndim; ++i)
         fX[i] = fMinX[i] + x[i] * fDeltaX[i];

      return (fFunc)(&fX[0]);
   }

   double  MinX(unsigned int i) { return fMinX[i]; }
   double  DeltaX(unsigned int i) { return fDeltaX[i]; }

private:

   const ROOT::Math::IMultiGenFunction & fFunc;
   std::vector<double> fX;
   std::vector<double> fMinX;
   std::vector<double> fDeltaX;

};



ClassImp(TFoamSampler)


//_______________________________________________________________________________
/**
   TFoamSampler class
   class implementing  the ROOT::Math::DistSampler interface using FOAM
   for sampling arbitrary distributions.


*/
TFoamSampler::TFoamSampler() : ROOT::Math::DistSampler(),
//    fOneDim(false)
//    fDiscrete(false),
//    fHasMode(false), fHasArea(false),
//    fMode(0), fArea(0),
   fFunc1D(0),
   fFoam(new TFoam("FOAM")  ),
   fFoamDist(0)
{}

TFoamSampler::~TFoamSampler() {
   assert(fFoam != 0);
   delete fFoam;
   if (fFoamDist) delete fFoamDist;
}

bool TFoamSampler::Init(const char *) {

   // initialize using default options
   ROOT::Math::DistSamplerOptions opt(0);
   ROOT::Math::IOptions * foamOpt  = ROOT::Math::DistSamplerOptions::FindDefault("Foam");
   if (foamOpt) opt.SetExtraOptions(*foamOpt);
   return Init(opt);
}

bool TFoamSampler::Init(const ROOT::Math::DistSamplerOptions & opt) {
   // initialize foam classes using the given algorithm
   assert (fFoam != 0 );
   if (NDim() == 0)  {
      Error("TFoamSampler::Init","Distribution function has not been set ! Need to call SetFunction first.");
      return false;
   }

   // initialize the foam
   fFoam->SetkDim(NDim() );

   // initialize random number
   if (!GetRandom()) SetRandom(gRandom);

   // create TFoamIntegrand class
   if (fFoamDist) delete fFoamDist;
   fFoamDist = new FoamDistribution(ParentPdf(),PdfRange());

   fFoam->SetRho(fFoamDist);
   // set print level
   fFoam->SetChat(opt.PrintLevel());

   // get extra options
   ROOT::Math::IOptions * fopt = opt.ExtraOptions();
   if (fopt) {
      int nval = 0;
      double fval = 0;
      if (fopt->GetIntValue("nCells", nval) ) fFoam->SetnCells(nval);
      if (fopt->GetIntValue("nCell1D", nval) && NDim() ==1) fFoam->SetnCells(nval);
      if (fopt->GetIntValue("nCellND", nval) && NDim()  >1) fFoam->SetnCells(nval);
      if (fopt->GetIntValue("nCell2D", nval) && NDim() ==2) fFoam->SetnCells(nval);
      if (fopt->GetIntValue("nCell3D", nval) && NDim() ==3) fFoam->SetnCells(nval);

      if (fopt->GetIntValue("nSample", nval) ) fFoam->SetnSampl(nval);
      if (fopt->GetIntValue("nBin", nval) ) fFoam->SetnBin(nval);
      if (fopt->GetIntValue("OptDrive",nval) ) fFoam->SetOptDrive(nval);
      if (fopt->GetIntValue("OptRej",nval) ) fFoam->SetOptRej(nval);
      if (fopt->GetRealValue("MaxWtRej",fval) ) fFoam->SetMaxWtRej(fval);


      if (fopt->GetIntValue("chatLevel", nval) ) fFoam->SetChat(nval);
   }
   fFoam->Initialize();

   return true;

}


void TFoamSampler::SetFunction(TF1 * pdf) {
   // set function from a TF1 pointer
   SetFunction<TF1>(*pdf, pdf->GetNdim());
}

void TFoamSampler::SetRandom(TRandom * r) {
   // set random generator (must be called before Init to have effect)
   fFoam->SetPseRan(r);
}

void TFoamSampler::SetSeed(unsigned int seed) {
   // set random generator seed (must be called before Init to have effect)
   TRandom * r = fFoam->GetPseRan();
   if (r) r->SetSeed(seed);
}

TRandom * TFoamSampler::GetRandom() {
   // get random generator used
   return  fFoam->GetPseRan();
}

// double TFoamSampler::Sample1D() {
//    // sample 1D distributions
//    return (fDiscrete) ? (double) fFoam->SampleDiscr() : fFoam->Sample();
// }

bool TFoamSampler::Sample(double * x) {
   // sample multi-dim distributions

   fFoam->MakeEvent();
   fFoam->GetMCvect(x);
   // adjust for the range
   for (unsigned int i = 0; i < NDim(); ++i)
      x[i] = ( (FoamDistribution*)fFoamDist)->MinX(i) + ( ( (FoamDistribution*) fFoamDist)->DeltaX(i))*x[i];

   return true;
}


bool TFoamSampler::SampleBin(double prob, double & value, double *error) {
   // sample a bin according to Poisson statistics

   TRandom * r = GetRandom();
   if (!r) return false;
   value = r->Poisson(prob);
   if (error) *error = std::sqrt(value);
   return true;
}
 TFoamSampler.cxx:1
 TFoamSampler.cxx:2
 TFoamSampler.cxx:3
 TFoamSampler.cxx:4
 TFoamSampler.cxx:5
 TFoamSampler.cxx:6
 TFoamSampler.cxx:7
 TFoamSampler.cxx:8
 TFoamSampler.cxx:9
 TFoamSampler.cxx:10
 TFoamSampler.cxx:11
 TFoamSampler.cxx:12
 TFoamSampler.cxx:13
 TFoamSampler.cxx:14
 TFoamSampler.cxx:15
 TFoamSampler.cxx:16
 TFoamSampler.cxx:17
 TFoamSampler.cxx:18
 TFoamSampler.cxx:19
 TFoamSampler.cxx:20
 TFoamSampler.cxx:21
 TFoamSampler.cxx:22
 TFoamSampler.cxx:23
 TFoamSampler.cxx:24
 TFoamSampler.cxx:25
 TFoamSampler.cxx:26
 TFoamSampler.cxx:27
 TFoamSampler.cxx:28
 TFoamSampler.cxx:29
 TFoamSampler.cxx:30
 TFoamSampler.cxx:31
 TFoamSampler.cxx:32
 TFoamSampler.cxx:33
 TFoamSampler.cxx:34
 TFoamSampler.cxx:35
 TFoamSampler.cxx:36
 TFoamSampler.cxx:37
 TFoamSampler.cxx:38
 TFoamSampler.cxx:39
 TFoamSampler.cxx:40
 TFoamSampler.cxx:41
 TFoamSampler.cxx:42
 TFoamSampler.cxx:43
 TFoamSampler.cxx:44
 TFoamSampler.cxx:45
 TFoamSampler.cxx:46
 TFoamSampler.cxx:47
 TFoamSampler.cxx:48
 TFoamSampler.cxx:49
 TFoamSampler.cxx:50
 TFoamSampler.cxx:51
 TFoamSampler.cxx:52
 TFoamSampler.cxx:53
 TFoamSampler.cxx:54
 TFoamSampler.cxx:55
 TFoamSampler.cxx:56
 TFoamSampler.cxx:57
 TFoamSampler.cxx:58
 TFoamSampler.cxx:59
 TFoamSampler.cxx:60
 TFoamSampler.cxx:61
 TFoamSampler.cxx:62
 TFoamSampler.cxx:63
 TFoamSampler.cxx:64
 TFoamSampler.cxx:65
 TFoamSampler.cxx:66
 TFoamSampler.cxx:67
 TFoamSampler.cxx:68
 TFoamSampler.cxx:69
 TFoamSampler.cxx:70
 TFoamSampler.cxx:71
 TFoamSampler.cxx:72
 TFoamSampler.cxx:73
 TFoamSampler.cxx:74
 TFoamSampler.cxx:75
 TFoamSampler.cxx:76
 TFoamSampler.cxx:77
 TFoamSampler.cxx:78
 TFoamSampler.cxx:79
 TFoamSampler.cxx:80
 TFoamSampler.cxx:81
 TFoamSampler.cxx:82
 TFoamSampler.cxx:83
 TFoamSampler.cxx:84
 TFoamSampler.cxx:85
 TFoamSampler.cxx:86
 TFoamSampler.cxx:87
 TFoamSampler.cxx:88
 TFoamSampler.cxx:89
 TFoamSampler.cxx:90
 TFoamSampler.cxx:91
 TFoamSampler.cxx:92
 TFoamSampler.cxx:93
 TFoamSampler.cxx:94
 TFoamSampler.cxx:95
 TFoamSampler.cxx:96
 TFoamSampler.cxx:97
 TFoamSampler.cxx:98
 TFoamSampler.cxx:99
 TFoamSampler.cxx:100
 TFoamSampler.cxx:101
 TFoamSampler.cxx:102
 TFoamSampler.cxx:103
 TFoamSampler.cxx:104
 TFoamSampler.cxx:105
 TFoamSampler.cxx:106
 TFoamSampler.cxx:107
 TFoamSampler.cxx:108
 TFoamSampler.cxx:109
 TFoamSampler.cxx:110
 TFoamSampler.cxx:111
 TFoamSampler.cxx:112
 TFoamSampler.cxx:113
 TFoamSampler.cxx:114
 TFoamSampler.cxx:115
 TFoamSampler.cxx:116
 TFoamSampler.cxx:117
 TFoamSampler.cxx:118
 TFoamSampler.cxx:119
 TFoamSampler.cxx:120
 TFoamSampler.cxx:121
 TFoamSampler.cxx:122
 TFoamSampler.cxx:123
 TFoamSampler.cxx:124
 TFoamSampler.cxx:125
 TFoamSampler.cxx:126
 TFoamSampler.cxx:127
 TFoamSampler.cxx:128
 TFoamSampler.cxx:129
 TFoamSampler.cxx:130
 TFoamSampler.cxx:131
 TFoamSampler.cxx:132
 TFoamSampler.cxx:133
 TFoamSampler.cxx:134
 TFoamSampler.cxx:135
 TFoamSampler.cxx:136
 TFoamSampler.cxx:137
 TFoamSampler.cxx:138
 TFoamSampler.cxx:139
 TFoamSampler.cxx:140
 TFoamSampler.cxx:141
 TFoamSampler.cxx:142
 TFoamSampler.cxx:143
 TFoamSampler.cxx:144
 TFoamSampler.cxx:145
 TFoamSampler.cxx:146
 TFoamSampler.cxx:147
 TFoamSampler.cxx:148
 TFoamSampler.cxx:149
 TFoamSampler.cxx:150
 TFoamSampler.cxx:151
 TFoamSampler.cxx:152
 TFoamSampler.cxx:153
 TFoamSampler.cxx:154
 TFoamSampler.cxx:155
 TFoamSampler.cxx:156
 TFoamSampler.cxx:157
 TFoamSampler.cxx:158
 TFoamSampler.cxx:159
 TFoamSampler.cxx:160
 TFoamSampler.cxx:161
 TFoamSampler.cxx:162
 TFoamSampler.cxx:163
 TFoamSampler.cxx:164
 TFoamSampler.cxx:165
 TFoamSampler.cxx:166
 TFoamSampler.cxx:167
 TFoamSampler.cxx:168
 TFoamSampler.cxx:169
 TFoamSampler.cxx:170
 TFoamSampler.cxx:171
 TFoamSampler.cxx:172
 TFoamSampler.cxx:173
 TFoamSampler.cxx:174
 TFoamSampler.cxx:175
 TFoamSampler.cxx:176
 TFoamSampler.cxx:177
 TFoamSampler.cxx:178
 TFoamSampler.cxx:179
 TFoamSampler.cxx:180
 TFoamSampler.cxx:181
 TFoamSampler.cxx:182
 TFoamSampler.cxx:183
 TFoamSampler.cxx:184
 TFoamSampler.cxx:185
 TFoamSampler.cxx:186
 TFoamSampler.cxx:187
 TFoamSampler.cxx:188
 TFoamSampler.cxx:189
 TFoamSampler.cxx:190
 TFoamSampler.cxx:191
 TFoamSampler.cxx:192
 TFoamSampler.cxx:193
 TFoamSampler.cxx:194
 TFoamSampler.cxx:195
 TFoamSampler.cxx:196
 TFoamSampler.cxx:197
 TFoamSampler.cxx:198
 TFoamSampler.cxx:199
 TFoamSampler.cxx:200
 TFoamSampler.cxx:201
 TFoamSampler.cxx:202
 TFoamSampler.cxx:203
 TFoamSampler.cxx:204
 TFoamSampler.cxx:205
 TFoamSampler.cxx:206
 TFoamSampler.cxx:207
 TFoamSampler.cxx:208