ROOT_528-00b_version: include/THnSparse.h Source File

00001 // @(#)root/hist:$Id: THnSparse.h 36894 2010-11-24 11:04:47Z moneta $
00002 // Author: Axel Naumann (2007-09-11)
00003 
00004 /*************************************************************************
00005  * Copyright (C) 1995-2007, Rene Brun and Fons Rademakers.               *
00006  * All rights reserved.                                                  *
00007  *                                                                       *
00008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
00009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
00010  *************************************************************************/
00011 
00012 #ifndef ROOT_THnSparse
00013 #define ROOT_THnSparse
00014 
00015 /*************************************************************************
00016 
00017  THnSparse: histogramming multi-dimensional, sparse distributions in
00018  a memory-efficient way.
00019 
00020 *************************************************************************/
00021 
00022 
00023 #ifndef ROOT_TExMap
00024 #include "TExMap.h"
00025 #endif
00026 #ifndef ROOT_TNamed
00027 #include "TNamed.h"
00028 #endif
00029 #ifndef ROOT_TObjArray
00030 #include "TObjArray.h"
00031 #endif
00032 #ifndef ROOT_TArrayD
00033 #include "TArrayD.h"
00034 #endif
00035 #ifndef ROOT_TFitResultPtr
00036 #include "TFitResultPtr.h"
00037 #endif
00038 
00039 // needed only for template instantiations of THnSparseT:
00040 #ifndef ROOT_TArrayF
00041 #include "TArrayF.h"
00042 #endif
00043 #ifndef ROOT_TArrayL
00044 #include "TArrayL.h"
00045 #endif
00046 #ifndef ROOT_TArrayI
00047 #include "TArrayI.h"
00048 #endif
00049 #ifndef ROOT_TArrayS
00050 #include "TArrayS.h"
00051 #endif
00052 #ifndef ROOT_TArrayC
00053 #include "TArrayC.h"
00054 #endif
00055 
00056 class TAxis;
00057 class TCollection;
00058 class TH1;
00059 class TH1D;
00060 class TH2D;
00061 class TH3D;
00062 class TF1;
00063 
00064 class THnSparseArrayChunk: public TObject {
00065  private:
00066 
00067    THnSparseArrayChunk(const THnSparseArrayChunk&); // Not implemented
00068    THnSparseArrayChunk& operator=(const THnSparseArrayChunk&); // Not implemented
00069 
00070  public:
00071    THnSparseArrayChunk():
00072       fCoordinateAllocationSize(-1), fSingleCoordinateSize(0), fCoordinatesSize(0), fCoordinates(0),
00073       fContent(0), fSumw2(0) {}
00074 
00075    THnSparseArrayChunk(Int_t coordsize, bool errors, TArray* cont);
00076    virtual ~THnSparseArrayChunk();
00077 
00078    Int_t    fCoordinateAllocationSize; //! size of the allocated coordinate buffer; -1 means none or fCoordinatesSize
00079    Int_t    fSingleCoordinateSize; // size of a single bin coordinate
00080    Int_t    fCoordinatesSize;      // size of the bin coordinate buffer
00081    Char_t  *fCoordinates;          //[fCoordinatesSize] compact bin coordinate buffer
00082    TArray  *fContent;              // bin content
00083    TArrayD *fSumw2;                // bin errors
00084 
00085    void AddBin(Int_t idx, const Char_t* idxbuf);
00086    void AddBinContent(Int_t idx, Double_t v = 1.) {
00087       fContent->SetAt(v + fContent->GetAt(idx), idx);
00088       if (fSumw2)
00089          fSumw2->SetAt(v * v+ fSumw2->GetAt(idx), idx);
00090    }
00091    void Sumw2();
00092    Int_t GetEntries() const { return fCoordinatesSize / fSingleCoordinateSize; }
00093    Bool_t Matches(Int_t idx, const Char_t* idxbuf) const {
00094       // Check whether bin at idx batches idxbuf.
00095       // If we don't store indexes we trust the caller that it does match,
00096       // see comment in THnSparseCompactBinCoord::GetHash().
00097       return fSingleCoordinateSize <= 8 ||
00098          !memcmp(fCoordinates + idx * fSingleCoordinateSize, idxbuf, fSingleCoordinateSize); }
00099 
00100    ClassDef(THnSparseArrayChunk, 1); // chunks of linearized bins
00101 };
00102 
00103 class THnSparseCompactBinCoord;
00104 
00105 class THnSparse: public TNamed {
00106  private:
00107    Int_t      fNdimensions;  // number of dimensions
00108    Int_t      fChunkSize;    // number of entries for each chunk
00109    Long64_t   fFilledBins;   // number of filled bins
00110    TObjArray  fAxes;         // axes of the histogram
00111    TObjArray  fBinContent;   // array of THnSparseArrayChunk
00112    TExMap     fBins;         //! filled bins
00113    TExMap     fBinsContinued;//! filled bins for non-unique hashes, containing pairs of (bin index 0, bin index 1)
00114    Double_t   fEntries;      // number of entries, spread over chunks
00115    Double_t   fTsumw;        // total sum of weights
00116    Double_t   fTsumw2;       // total sum of weights squared; -1 if no errors are calculated
00117    TArrayD    fTsumwx;       // total sum of weight*X for each dimension
00118    TArrayD    fTsumwx2;      // total sum of weight*X*X for each dimension
00119    THnSparseCompactBinCoord *fCompactCoord; //! compact coordinate
00120    Double_t  *fIntegral;     //! array with bin weight sums
00121    enum {
00122       kNoInt,
00123       kValidInt,
00124       kInvalidInt
00125    } fIntegralStatus;        //! status of integral
00126 
00127    THnSparse(const THnSparse&); // Not implemented
00128    THnSparse& operator=(const THnSparse&); // Not implemented
00129 
00130  protected:
00131 
00132    THnSparse();
00133    THnSparse(const char* name, const char* title, Int_t dim,
00134              const Int_t* nbins, const Double_t* xmin, const Double_t* xmax,
00135              Int_t chunksize);
00136    Int_t GetChunkSize() const { return fChunkSize; }
00137    THnSparseCompactBinCoord* GetCompactCoord() const;
00138    THnSparseArrayChunk* GetChunk(Int_t idx) const {
00139       return (THnSparseArrayChunk*) fBinContent[idx]; }
00140 
00141    THnSparseArrayChunk* AddChunk();
00142    void FillExMap();
00143    virtual TArray* GenerateArray() const = 0;
00144    Long64_t GetBinIndexForCurrentBin(Bool_t allocate);
00145    Long64_t Fill(Long64_t bin, Double_t w) {
00146       // Increment the bin content of "bin" by "w",
00147       // return the bin index.
00148       fEntries += 1;
00149       if (GetCalculateErrors()) {
00150          fTsumw += w;
00151          fTsumw2 += w*w;
00152       }
00153       fIntegralStatus = kInvalidInt;
00154       THnSparseArrayChunk* chunk = GetChunk(bin / fChunkSize);
00155       chunk->AddBinContent(bin % fChunkSize, w);
00156       return bin;
00157    }
00158    THnSparse* CloneEmpty(const char* name, const char* title,
00159                          const TObjArray* axes, Int_t chunksize,
00160                          Bool_t keepTargetAxis) const;
00161 
00162    Bool_t CheckConsistency(const THnSparse *h, const char *tag) const;
00163    Bool_t IsInRange(Int_t *coord) const;
00164    TH1* CreateHist(const char* name, const char* title,
00165                    const TObjArray* axes, Bool_t keepTargetAxis) const;
00166    TObject* ProjectionAny(Int_t ndim, const Int_t* dim,
00167                           Bool_t wantSparse, Option_t* option = "") const;
00168    Bool_t PrintBin(Long64_t idx, Int_t* coord, Option_t* options) const;
00169    void AddInternal(const THnSparse* h, Double_t c, Bool_t rebinned);
00170 
00171  public:
00172    virtual ~THnSparse();
00173 
00174    static THnSparse* CreateSparse(const char* name, const char* title,
00175                                   const TH1* h1, Int_t ChunkSize = 1024 * 16);
00176 
00177    Int_t  GetNChunks() const { return fBinContent.GetEntriesFast(); }
00178    TObjArray* GetListOfAxes() { return &fAxes; }
00179    TAxis* GetAxis(Int_t dim) const { return (TAxis*)fAxes[dim]; }
00180 
00181    Long64_t Fill(const Double_t *x, Double_t w = 1.) {
00182       if (GetCalculateErrors()) {
00183          for (Int_t d = 0; d < fNdimensions; ++d) {
00184             const Double_t xd = x[d];
00185             fTsumwx[d]  += w * xd;
00186             fTsumwx2[d] += w * xd * xd;
00187          }
00188       }
00189       return Fill(GetBin(x), w);
00190    }
00191    Long64_t Fill(const char* name[], Double_t w = 1.) {
00192       return Fill(GetBin(name), w);
00193    }
00194 
00195    TFitResultPtr Fit(TF1 *f1 ,Option_t *option = "", Option_t *goption = "");
00196    TList* GetListOfFunctions() { return 0; }
00197 
00198    Double_t GetEntries() const { return fEntries; }
00199    Double_t GetWeightSum() const { return fTsumw; }
00200    Long64_t GetNbins() const { return fFilledBins; }
00201    Int_t    GetNdimensions() const { return fNdimensions; }
00202    Bool_t   GetCalculateErrors() const { return fTsumw2 >= 0.; }
00203    void     CalculateErrors(Bool_t calc = kTRUE) {
00204       // Calculate errors (or not if "calc" == kFALSE)
00205       if (calc) Sumw2();
00206       else fTsumw2 = -1.;
00207    }
00208 
00209    Long64_t GetBin(const Int_t* idx, Bool_t allocate = kTRUE);
00210    Long64_t GetBin(const Double_t* x, Bool_t allocate = kTRUE);
00211    Long64_t GetBin(const char* name[], Bool_t allocate = kTRUE);
00212 
00213    void SetBinEdges(Int_t idim, const Double_t* bins);
00214    void SetBinContent(const Int_t* x, Double_t v);
00215    void SetBinContent(Long64_t bin, Double_t v);
00216    void SetBinError(const Int_t* x, Double_t e);
00217    void SetBinError(Long64_t bin, Double_t e);
00218    void AddBinContent(const Int_t* x, Double_t v = 1.);
00219    void AddBinContent(Long64_t bin, Double_t v = 1.);
00220    void SetEntries(Double_t entries) { fEntries = entries; }
00221    void SetTitle(const char *title);
00222 
00223    Double_t GetBinContent(const Int_t *idx) const;
00224    Double_t GetBinContent(Long64_t bin, Int_t* idx = 0) const;
00225    Double_t GetBinError(const Int_t *idx) const;
00226    Double_t GetBinError(Long64_t linidx) const;
00227 
00228    Double_t GetSparseFractionBins() const;
00229    Double_t GetSparseFractionMem() const;
00230 
00231    Double_t GetSumw() const  { return fTsumw; }
00232    Double_t GetSumw2() const { return fTsumw2; }
00233    Double_t GetSumwx(Int_t dim) const  { return fTsumwx[dim]; }
00234    Double_t GetSumwx2(Int_t dim) const { return fTsumwx2[dim]; }
00235 
00236    TH1D*      Projection(Int_t xDim, Option_t* option = "") const;
00237    TH2D*      Projection(Int_t yDim, Int_t xDim,
00238                          Option_t* option = "") const;
00239    TH3D*      Projection(Int_t xDim, Int_t yDim, Int_t zDim,
00240                          Option_t* option = "") const;
00241    THnSparse* Projection(Int_t ndim, const Int_t* dim,
00242                          Option_t* option = "") const;
00243 
00244    THnSparse* Rebin(Int_t group) const;
00245    THnSparse* Rebin(const Int_t* group) const;
00246 
00247    Long64_t   Merge(TCollection* list);
00248 
00249    void Scale(Double_t c);
00250    void Add(const THnSparse* h, Double_t c=1.);
00251    void Multiply(const THnSparse* h);
00252    void Multiply(TF1* f, Double_t c = 1.);
00253    void Divide(const THnSparse* h);
00254    void Divide(const THnSparse* h1, const THnSparse* h2, Double_t c1 = 1., Double_t c2 = 1., Option_t* option="");
00255    void RebinnedAdd(const THnSparse* h, Double_t c=1.);
00256 
00257    void Reset(Option_t* option = "");
00258    void Sumw2();
00259 
00260    Double_t ComputeIntegral();
00261    void GetRandom(Double_t *rand, Bool_t subBinRandom = kTRUE);
00262 
00263    void Print(Option_t* option = "") const;
00264    void PrintEntries(Long64_t from = 0, Long64_t howmany = -1, Option_t* options = 0) const;
00265    void PrintBin(Int_t* coord, Option_t* options) const {
00266       PrintBin(-1, coord, options);
00267    }
00268    void PrintBin(Long64_t idx, Option_t* options) const;
00269 
00270    //void Draw(Option_t* option = "");
00271 
00272    ClassDef(THnSparse, 2); // Interfaces of sparse n-dimensional histogram
00273 };
00274 
00275 
00276 
00277 //______________________________________________________________________________
00278 //
00279 // Templated implementation of the abstract base THnSparse.
00280 // All functionality and the interfaces to be used are in THnSparse!
00281 //
00282 // THnSparse does not know how to store any bin content itself. Instead, this
00283 // is delegated to the derived, templated class: the template parameter decides
00284 // what the format for the bin content is. In fact it even defines the array
00285 // itself; possible implementations probably derive from TArray.
00286 //
00287 // Typedefs exist for template parematers with ROOT's generic types:
00288 //
00289 //   Templated name        Typedef       Bin content type
00290 //   THnSparseT<TArrayC>   THnSparseC    Char_r
00291 //   THnSparseT<TArrayS>   THnSparseS    Short_t
00292 //   THnSparseT<TArrayI>   THnSparseI    Int_t
00293 //   THnSparseT<TArrayL>   THnSparseL    Long_t
00294 //   THnSparseT<TArrayF>   THnSparseF    Float_t
00295 //   THnSparseT<TArrayD>   THnSparseD    Double_t
00296 //
00297 // We recommend to use THnSparseC wherever possible, and to map its value space
00298 // of 256 possible values to e.g. float values outside the class. This saves an
00299 // enourmous amount of memory. Only if more than 256 values need to be
00300 // distinguished should e.g. THnSparseS or even THnSparseF be chosen.
00301 //
00302 // Implementation detail: the derived, templated class is kept extremely small
00303 // on purpose. That way the (templated thus inlined) uses of this class will
00304 // only create a small amount of machine code, in contrast to e.g. STL.
00305 //______________________________________________________________________________
00306 
00307 template <class CONT>
00308 class THnSparseT: public THnSparse {
00309  public:
00310    THnSparseT() {}
00311    THnSparseT(const char* name, const char* title, Int_t dim,
00312               const Int_t* nbins, const Double_t* xmin = 0,
00313               const Double_t* xmax = 0, Int_t chunksize = 1024 * 16):
00314       THnSparse(name, title, dim, nbins, xmin, xmax, chunksize) {}
00315 
00316    TArray* GenerateArray() const { return new CONT(GetChunkSize()); }
00317  private:
00318    ClassDef(THnSparseT, 1); // Sparse n-dimensional histogram with templated content
00319 };
00320 
00321 typedef THnSparseT<TArrayD> THnSparseD;
00322 typedef THnSparseT<TArrayF> THnSparseF;
00323 typedef THnSparseT<TArrayL> THnSparseL;
00324 typedef THnSparseT<TArrayI> THnSparseI;
00325 typedef THnSparseT<TArrayS> THnSparseS;
00326 typedef THnSparseT<TArrayC> THnSparseC;
00327 
00328 #endif //  ROOT_THnSparse