Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
THnSparse.cxx
Go to the documentation of this file.
1// @(#)root/hist:$Id$
2// Author: Axel Naumann (2007-09-11)
3
4/*************************************************************************
5 * Copyright (C) 1995-2012, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12#include "THnSparse.h"
13
14#include "TAxis.h"
15#include "TClass.h"
16#include "TDataMember.h"
17#include "TDataType.h"
18
19namespace {
20//______________________________________________________________________________
21//
22// THnSparseBinIter iterates over all filled bins of a THnSparse.
23//______________________________________________________________________________
24
25 class THnSparseBinIter: public ROOT::Internal::THnBaseBinIter {
26 public:
27 THnSparseBinIter(Bool_t respectAxisRange, const THnSparse* hist):
28 ROOT::Internal::THnBaseBinIter(respectAxisRange), fHist(hist),
29 fNbins(hist->GetNbins()), fIndex(-1) {
30 // Construct a THnSparseBinIter
31 fCoord = new Int_t[hist->GetNdimensions()];
32 fCoord[0] = -1;
33 }
34 virtual ~THnSparseBinIter() { delete [] fCoord; }
35
36 virtual Int_t GetCoord(Int_t dim) const;
37 virtual Long64_t Next(Int_t* coord = 0);
38
39 private:
40 THnSparseBinIter(const THnSparseBinIter&); // intentionally unimplemented
41 THnSparseBinIter& operator=(const THnSparseBinIter&); // intentionally unimplemented
42
43 const THnSparse* fHist;
44 Int_t* fCoord; // coord buffer for fIndex; fCoord[0] == -1 if not yet calculated
45 Long64_t fNbins; // number of bins to iterate over
46 Long64_t fIndex; // current bin index
47 };
48}
49
50Int_t THnSparseBinIter::GetCoord(Int_t dim) const
51{
52 if (fCoord[0] == -1) {
53 fHist->GetBinContent(fIndex, fCoord);
54 }
55 return fCoord[dim];
56}
57
58Long64_t THnSparseBinIter::Next(Int_t* coord /*= 0*/)
59{
60 // Get next bin index (in range if RespectsAxisRange()).
61 // If coords != 0, set it to the index's axis coordinates
62 // (i.e. coord must point to an array of Int_t[fNdimension]
63 if (!fHist) return -1;
64
65 fCoord[0] = -1;
66 Int_t* useCoordBuf = fCoord;
67 if (coord) {
68 useCoordBuf = coord;
69 coord[0] = -1;
70 }
71
72 do {
73 ++fIndex;
74 if (fIndex >= fHist->GetNbins()) {
75 fHist = 0;
76 return -1;
77 }
78 if (RespectsAxisRange()) {
79 fHist->GetBinContent(fIndex, useCoordBuf);
80 }
81 } while (RespectsAxisRange()
82 && !fHist->IsInRange(useCoordBuf)
83 && (fHaveSkippedBin = kTRUE /* assignment! */));
84
85 if (coord && coord[0] == -1) {
86 if (fCoord[0] == -1) {
87 fHist->GetBinContent(fIndex, coord);
88 } else {
89 memcpy(coord, fCoord, fHist->GetNdimensions() * sizeof(Int_t));
90 }
91 }
92
93 return fIndex;
94}
95
96
97
98/** \class THnSparseCoordCompression
99THnSparseCoordCompression is a class used by THnSparse internally. It
100represents a compacted n-dimensional array of bin coordinates (indices).
101As the total number of bins in each dimension is known by THnSparse, bin
102indices can be compacted to only use the amount of bins needed by the total
103number of bins in each dimension. E.g. for a THnSparse with
104{15, 100, 2, 20, 10, 100} bins per dimension, a bin index will only occupy
10528 bits (4+7+1+5+4+7), i.e. less than a 32bit integer. The tricky part is
106the fast compression and decompression, the platform-independent storage
107(think of endianness: the bits of the number 0x123456 depend on the
108platform), and the hashing needed by THnSparseArrayChunk.
109*/
110
111
112class THnSparseCoordCompression {
113public:
114 THnSparseCoordCompression(Int_t dim, const Int_t* nbins);
115 THnSparseCoordCompression(const THnSparseCoordCompression& other);
116 ~THnSparseCoordCompression();
117
118 THnSparseCoordCompression& operator=(const THnSparseCoordCompression& other);
119
120 ULong64_t GetHashFromBuffer(const Char_t* buf) const;
121 Int_t GetBufferSize() const { return fCoordBufferSize; }
122 Int_t GetNdimensions() const { return fNdimensions; }
123 void SetCoordFromBuffer(const Char_t* buf_in, Int_t* coord_out) const;
124 ULong64_t SetBufferFromCoord(const Int_t* coord_in, Char_t* buf_out) const;
125
126protected:
127 Int_t GetNumBits(Int_t n) const {
128 // return the number of bits allocated by the number "n"
129 Int_t r = (n > 0);
130 while (n/=2) ++r;
131 return r;
132 }
133private:
134 Int_t fNdimensions; // number of dimensions
135 Int_t fCoordBufferSize; // size of coordbuf
136 Int_t *fBitOffsets; //[fNdimensions + 1] bit offset of each axis index
137};
138
139
140//______________________________________________________________________________
141//______________________________________________________________________________
142
143
144////////////////////////////////////////////////////////////////////////////////
145/// Initialize a THnSparseCoordCompression object with "dim" dimensions
146/// and "bins" holding the number of bins for each dimension; it
147/// stores the
148
149THnSparseCoordCompression::THnSparseCoordCompression(Int_t dim, const Int_t* nbins):
150 fNdimensions(dim), fCoordBufferSize(0), fBitOffsets(0)
151{
152 fBitOffsets = new Int_t[dim + 1];
153
154 int shift = 0;
155 for (Int_t i = 0; i < dim; ++i) {
156 fBitOffsets[i] = shift;
157 shift += GetNumBits(nbins[i] + 2);
158 }
159 fBitOffsets[dim] = shift;
160 fCoordBufferSize = (shift + 7) / 8;
161}
162
163
164////////////////////////////////////////////////////////////////////////////////
165/// Construct a THnSparseCoordCompression from another one
166
167THnSparseCoordCompression::THnSparseCoordCompression(const THnSparseCoordCompression& other)
168{
169 fNdimensions = other.fNdimensions;
170 fCoordBufferSize = other.fCoordBufferSize;
171 fBitOffsets = new Int_t[fNdimensions + 1];
172 memcpy(fBitOffsets, other.fBitOffsets, sizeof(Int_t) * fNdimensions);
173}
174
175
176////////////////////////////////////////////////////////////////////////////////
177/// Set this to other if different.
178
179THnSparseCoordCompression& THnSparseCoordCompression::operator=(const THnSparseCoordCompression& other)
180{
181 if (&other == this) return *this;
182
183 fNdimensions = other.fNdimensions;
184 fCoordBufferSize = other.fCoordBufferSize;
185 delete [] fBitOffsets;
186 fBitOffsets = new Int_t[fNdimensions + 1];
187 memcpy(fBitOffsets, other.fBitOffsets, sizeof(Int_t) * fNdimensions);
188 return *this;
189}
190
191
192////////////////////////////////////////////////////////////////////////////////
193/// destruct a THnSparseCoordCompression
194
195THnSparseCoordCompression::~THnSparseCoordCompression()
196{
197 delete [] fBitOffsets;
198}
199
200
201////////////////////////////////////////////////////////////////////////////////
202/// Given the compressed coordinate buffer buf_in, calculate ("decompact")
203/// the bin coordinates and return them in coord_out.
204
205void THnSparseCoordCompression::SetCoordFromBuffer(const Char_t* buf_in,
206 Int_t* coord_out) const
207{
208 for (Int_t i = 0; i < fNdimensions; ++i) {
209 const Int_t offset = fBitOffsets[i] / 8;
210 Int_t shift = fBitOffsets[i] % 8;
211 Int_t nbits = fBitOffsets[i + 1] - fBitOffsets[i];
212 const UChar_t* pbuf = (const UChar_t*) buf_in + offset;
213 coord_out[i] = *pbuf >> shift;
214 Int_t subst = (Int_t) -1;
215 subst = subst << nbits;
216 nbits -= (8 - shift);
217 shift = 8 - shift;
218 for (Int_t n = 0; n * 8 < nbits; ++n) {
219 ++pbuf;
220 coord_out[i] += *pbuf << shift;
221 shift += 8;
222 }
223 coord_out[i] &= ~subst;
224 }
225}
226
227
228////////////////////////////////////////////////////////////////////////////////
229/// Given the cbin coordinates coord_in, calculate ("compact")
230/// the bin coordinates and return them in buf_in.
231/// Return the hash value.
232
233ULong64_t THnSparseCoordCompression::SetBufferFromCoord(const Int_t* coord_in,
234 Char_t* buf_out) const
235{
236 if (fCoordBufferSize <= 8) {
237 ULong64_t l64buf = 0;
238 for (Int_t i = 0; i < fNdimensions; ++i) {
239 l64buf += ((ULong64_t)((UInt_t)coord_in[i])) << fBitOffsets[i];
240 }
241 memcpy(buf_out, &l64buf, sizeof(Long64_t));
242 return l64buf;
243 }
244
245 // else: doesn't fit into a Long64_t:
246 memset(buf_out, 0, fCoordBufferSize);
247 for (Int_t i = 0; i < fNdimensions; ++i) {
248 const Int_t offset = fBitOffsets[i] / 8;
249 const Int_t shift = fBitOffsets[i] % 8;
250 ULong64_t val = coord_in[i];
251
252 Char_t* pbuf = buf_out + offset;
253 *pbuf += 0xff & (val << shift);
254 val = val >> (8 - shift);
255 while (val) {
256 ++pbuf;
257 *pbuf += 0xff & val;
258 val = val >> 8;
259 }
260 }
261
262 return GetHashFromBuffer(buf_out);
263}
264
265/*
266////////////////////////////////////////////////////////////////////////////////
267/// Calculate hash from bin indexes.
268
269ULong64_t THnSparseCoordCompression::GetHashFromCoords(const Int_t* coord) const
270{
271 // Bins are addressed in two different modes, depending
272 // on whether the compact bin index fits into a Long64_t or not.
273 // If it does, we can use it as a "perfect hash" for the TExMap.
274 // If not we build a hash from the compact bin index, and use that
275 // as the TExMap's hash.
276
277 if (fCoordBufferSize <= 8) {
278 // fits into a Long64_t
279 ULong64_t hash1 = 0;
280 for (Int_t i = 0; i < fNdimensions; ++i) {
281 hash1 += coord[i] << fBitOffsets[i];
282 }
283 return hash1;
284 }
285
286 // else: doesn't fit into a Long64_t:
287 memset(coord, 0, fCoordBufferSize);
288 for (Int_t i = 0; i < fNdimensions; ++i) {
289 const Int_t offset = fBitOffsets[i] / 8;
290 const Int_t shift = fBitOffsets[i] % 8;
291 ULong64_t val = coord[i];
292
293 Char_t* pbuf = fCoordBuffer + offset;
294 *pbuf += 0xff & (val << shift);
295 val = val >> (8 - shift);
296 while (val) {
297 ++pbuf;
298 *pbuf += 0xff & val;
299 val = val >> 8;
300 }
301 }
302
303 ULong64_t hash = 5381;
304 Char_t* str = fCoordBuffer;
305 while (str - fCoordBuffer < fCoordBufferSize) {
306 hash *= 5;
307 hash += *(str++);
308 }
309 return hash;
310}
311*/
312
313
314////////////////////////////////////////////////////////////////////////////////
315/// Calculate hash from compact bin index.
316
317ULong64_t THnSparseCoordCompression::GetHashFromBuffer(const Char_t* buf) const
318{
319 // Bins are addressed in two different modes, depending
320 // on whether the compact bin index fits into a Long64_t or not.
321 // If it does, we can use it as a "perfect hash" for the TExMap.
322 // If not we build a hash from the compact bin index, and use that
323 // as the TExMap's hash.
324
325 if (fCoordBufferSize <= 8) {
326 // fits into a Long64_t
327 ULong64_t hash1 = 0;
328 memcpy(&hash1, buf, fCoordBufferSize);
329 return hash1;
330 }
331
332 // else: doesn't fit into a Long64_t:
333 ULong64_t hash = 5381;
334 const Char_t* str = buf;
335 while (str - buf < fCoordBufferSize) {
336 hash *= 5;
337 hash += *(str++);
338 }
339 return hash;
340}
341
342
343
344
345/** \class THnSparseCompactBinCoord
346THnSparseCompactBinCoord is a class used by THnSparse internally. It
347maps between an n-dimensional array of bin coordinates (indices) and
348its compact version, the THnSparseCoordCompression.
349*/
350
351class THnSparseCompactBinCoord: public THnSparseCoordCompression {
352public:
353 THnSparseCompactBinCoord(Int_t dim, const Int_t* nbins);
354 ~THnSparseCompactBinCoord();
355 Int_t* GetCoord() { return fCurrentBin; }
356 const Char_t* GetBuffer() const { return fCoordBuffer; }
357 ULong64_t GetHash() const { return fHash; }
358 void UpdateCoord() {
359 fHash = SetBufferFromCoord(fCurrentBin, fCoordBuffer);
360 }
361 void SetCoord(const Int_t* coord) {
362 memcpy(fCurrentBin, coord, sizeof(Int_t) * GetNdimensions());
363 fHash = SetBufferFromCoord(coord, fCoordBuffer);
364 }
365 void SetBuffer(const Char_t* buf) {
366 memcpy(fCoordBuffer, buf, GetBufferSize());
367 fHash = GetHashFromBuffer(fCoordBuffer);
368 }
369
370private:
371 // intentionally not implemented
372 THnSparseCompactBinCoord(const THnSparseCompactBinCoord&);
373 // intentionally not implemented
374 THnSparseCompactBinCoord& operator=(const THnSparseCompactBinCoord&);
375
376private:
377 ULong64_t fHash; // hash for current coordinates; 0 if not calculated
378 Char_t *fCoordBuffer; // compact buffer of coordinates
379 Int_t *fCurrentBin; // current coordinates
380};
381
382
383//______________________________________________________________________________
384//______________________________________________________________________________
385
386
387////////////////////////////////////////////////////////////////////////////////
388/// Initialize a THnSparseCompactBinCoord object with "dim" dimensions
389/// and "bins" holding the number of bins for each dimension.
390
391THnSparseCompactBinCoord::THnSparseCompactBinCoord(Int_t dim, const Int_t* nbins):
392 THnSparseCoordCompression(dim, nbins),
393 fHash(0), fCoordBuffer(0), fCurrentBin(0)
394{
395 fCurrentBin = new Int_t[dim];
396 size_t bufAllocSize = GetBufferSize();
397 if (bufAllocSize < sizeof(Long64_t))
398 bufAllocSize = sizeof(Long64_t);
399 fCoordBuffer = new Char_t[bufAllocSize];
400}
401
402
403////////////////////////////////////////////////////////////////////////////////
404/// destruct a THnSparseCompactBinCoord
405
406THnSparseCompactBinCoord::~THnSparseCompactBinCoord()
407{
408 delete [] fCoordBuffer;
409 delete [] fCurrentBin;
410}
411
412/** \class THnSparseArrayChunk
413THnSparseArrayChunk is used internally by THnSparse.
414THnSparse stores its (dynamic size) array of bin coordinates and their
415contents (and possibly errors) in a TObjArray of THnSparseArrayChunk. Each
416of the chunks holds an array of THnSparseCompactBinCoord and the content
417(a TArray*), which is created outside (by the templated derived classes of
418THnSparse) and passed in at construction time.
419*/
420
422
423////////////////////////////////////////////////////////////////////////////////
424/// (Default) initialize a chunk. Takes ownership of cont (~THnSparseArrayChunk deletes it),
425/// and create an ArrayF for errors if "errors" is true.
426
428 fCoordinateAllocationSize(-1), fSingleCoordinateSize(coordsize), fCoordinatesSize(0),
429 fCoordinates(0), fContent(cont),
430 fSumw2(0)
431{
434 if (errors) Sumw2();
435}
436
437////////////////////////////////////////////////////////////////////////////////
438/// Destructor
439
441{
442 delete fContent;
443 delete [] fCoordinates;
444 delete fSumw2;
445}
446
447////////////////////////////////////////////////////////////////////////////////
448/// Create a new bin in this chunk
449
450void THnSparseArrayChunk::AddBin(Int_t idx, const Char_t* coordbuf)
451{
452 // When streaming out only the filled chunk is saved.
453 // When reading back only the memory needed for that filled part gets
454 // allocated. We need to check whether the allowed chunk size is
455 // bigger than the allocated size. If fCoordinateAllocationSize is
456 // set to -1 this chunk has been allocated by the streamer and the
457 // buffer allocation size is defined by [fCoordinatesSize]. In that
458 // case we need to compare fCoordinatesSize to
459 // fSingleCoordinateSize * fContent->GetSize()
460 // to determine whether we need to expand the buffer.
461 if (fCoordinateAllocationSize == -1 && fContent) {
463 if (fCoordinatesSize < chunksize) {
464 // need to re-allocate:
465 Char_t *newcoord = new Char_t[chunksize];
466 memcpy(newcoord, fCoordinates, fCoordinatesSize);
467 delete [] fCoordinates;
468 fCoordinates = newcoord;
469 }
470 fCoordinateAllocationSize = chunksize;
471 }
472
475}
476
477////////////////////////////////////////////////////////////////////////////////
478/// Turn on support of errors
479
481{
482 if (!fSumw2)
483 fSumw2 = new TArrayD(fContent->GetSize());
484 // fill the structure with the current content
485 for (Int_t bin=0; bin < fContent->GetSize(); bin++) {
486 fSumw2->fArray[bin] = fContent->GetAt(bin);
487 }
488
489}
490
491
492/** \class THnSparse
493 \ingroup Hist
494
495Efficient multidimensional histogram.
496
497Use a THnSparse instead of TH1 / TH2 / TH3 / array for histogramming when
498only a small fraction of bins is filled. A 10-dimensional histogram with 10
499bins per dimension has 10^10 bins; in a naive implementation this will not
500fit in memory. THnSparse only allocates memory for the bins that have
501non-zero bin content instead, drastically reducing both the memory usage
502and the access time.
503
504To construct a THnSparse object you must use one of its templated, derived
505classes:
506- THnSparseD (typedef for THnSparseT<ArrayD>): bin content held by a Double_t,
507- THnSparseF (typedef for THnSparseT<ArrayF>): bin content held by a Float_t,
508- THnSparseL (typedef for THnSparseT<ArrayL>): bin content held by a Long_t,
509- THnSparseI (typedef for THnSparseT<ArrayI>): bin content held by an Int_t,
510- THnSparseS (typedef for THnSparseT<ArrayS>): bin content held by a Short_t,
511- THnSparseC (typedef for THnSparseT<ArrayC>): bin content held by a Char_t,
512
513They take name and title, the number of dimensions, and for each dimension
514the number of bins, the minimal, and the maximal value on the dimension's
515axis. A TH2 h("h","h",10, 0., 10., 20, -5., 5.) would correspond to
516
517 Int_t bins[2] = {10, 20};
518 Double_t xmin[2] = {0., -5.};
519 Double_t xmax[2] = {10., 5.};
520 THnSparseD hs("hs", "hs", 2, bins, xmin, xmax);
521
522## Filling
523A THnSparse is filled just like a regular histogram, using
524THnSparse::Fill(x, weight), where x is a n-dimensional Double_t value.
525To take errors into account, Sumw2() must be called before filling the
526histogram.
527
528Bins are allocated as needed; the status of the allocation can be observed
529by GetSparseFractionBins(), GetSparseFractionMem().
530
531## Fast Bin Content Access
532When iterating over a THnSparse one should only look at filled bins to save
533processing time. The number of filled bins is returned by
534THnSparse::GetNbins(); the bin content for each (linear) bin number can
535be retrieved by THnSparse::GetBinContent(linidx, (Int_t*)coord).
536After the call, coord will contain the bin coordinate of each axis for the bin
537with linear index linidx. A possible call would be
538
539 std::cout << hs.GetBinContent(0, coord);
540 std::cout <<" is the content of bin [x = " << coord[0] "
541 << " | y = " << coord[1] << "]" << std::endl;
542
543## Efficiency
544TH1 and TH2 are generally faster than THnSparse for one and two dimensional
545distributions. THnSparse becomes competitive for a sparsely filled TH3
546with large numbers of bins per dimension. The tutorial sparsehist.C
547shows the turning point. On a AMD64 with 8GB memory, THnSparse "wins"
548starting with a TH3 with 30 bins per dimension. Using a THnSparse for a
549one-dimensional histogram is only reasonable if it has a huge number of bins.
550
551## Projections
552The dimensionality of a THnSparse can be reduced by projecting it to
5531, 2, 3, or n dimensions, which can be represented by a TH1, TH2, TH3, or
554a THnSparse. See the Projection() members. To only project parts of the
555histogram, call
556
557 THnSparse::GetAxis(12)->SetRange(from_bin, to_bin);
558
559## Internal Representation
560An entry for a filled bin consists of its n-dimensional coordinates and
561its bin content. The coordinates are compacted to use as few bits as
562possible; e.g. a histogram with 10 bins in x and 20 bins in y will only
563use 4 bits for the x representation and 5 bits for the y representation.
564This is handled by the internal class THnSparseCompactBinCoord.
565Bin data (content and coordinates) are allocated in chunks of size
566fChunkSize; this parameter can be set when constructing a THnSparse. Each
567chunk is represented by an object of class THnSparseArrayChunk.
568
569Translation from an n-dimensional bin coordinate to the linear index within
570the chunks is done by GetBin(). It creates a hash from the compacted bin
571coordinates (the hash of a bin coordinate is the compacted coordinate itself
572if it takes less than 8 bytes, the size of a Long64_t.
573This hash is used to lookup the linear index in the TExMap member fBins;
574the coordinates of the entry fBins points to is compared to the coordinates
575passed to GetBin(). If they do not match, these two coordinates have the same
576hash - which is extremely unlikely but (for the case where the compact bin
577coordinates are larger than 4 bytes) possible. In this case, fBinsContinued
578contains a chain of linear indexes with the same hash. Iterating through this
579chain and comparing each bin coordinates with the one passed to GetBin() will
580retrieve the matching bin.
581*/
582
583
585
586////////////////////////////////////////////////////////////////////////////////
587/// Construct an empty THnSparse.
588
590 fChunkSize(1024), fFilledBins(0), fCompactCoord(0)
591{
593}
594
595////////////////////////////////////////////////////////////////////////////////
596/// Construct a THnSparse with "dim" dimensions,
597/// with chunksize as the size of the chunks.
598/// "nbins" holds the number of bins for each dimension;
599/// "xmin" and "xmax" the minimal and maximal value for each dimension.
600/// The arrays "xmin" and "xmax" can be NULL; in that case SetBinEdges()
601/// must be called for each dimension.
602
603THnSparse::THnSparse(const char* name, const char* title, Int_t dim,
604 const Int_t* nbins, const Double_t* xmin, const Double_t* xmax,
605 Int_t chunksize):
606 THnBase(name, title, dim, nbins, xmin, xmax),
607 fChunkSize(chunksize), fFilledBins(0), fCompactCoord(0)
608{
609 fCompactCoord = new THnSparseCompactBinCoord(dim, nbins);
611}
612
613////////////////////////////////////////////////////////////////////////////////
614/// Destruct a THnSparse
615
617 delete fCompactCoord;
618}
619
620////////////////////////////////////////////////////////////////////////////////
621/// Add "v" to the content of bin with index "bin"
622
624{
626 bin %= fChunkSize;
627 v += chunk->fContent->GetAt(bin);
628 return chunk->fContent->SetAt(v, bin);
629}
630
631////////////////////////////////////////////////////////////////////////////////
632/// Create a new chunk of bin content
633
635{
636 THnSparseArrayChunk* chunk =
637 new THnSparseArrayChunk(GetCompactCoord()->GetBufferSize(),
639 fBinContent.AddLast(chunk);
640 return chunk;
641}
642
643////////////////////////////////////////////////////////////////////////////////
644/// Initialize the storage of a histogram created via Init()
645
646void THnSparse::InitStorage(Int_t* nbins, Int_t chunkSize)
647{
648 fChunkSize = chunkSize;
649 fCompactCoord = new THnSparseCompactBinCoord(fNdimensions, nbins);
650}
651
652////////////////////////////////////////////////////////////////////////////////
653///We have been streamed; set up fBins
654
656{
657 TIter iChunk(&fBinContent);
658 THnSparseArrayChunk* chunk = 0;
659 THnSparseCoordCompression compactCoord(*GetCompactCoord());
660 Long64_t idx = 0;
661 if (2 * GetNbins() > fBins.Capacity())
662 fBins.Expand(3 * GetNbins());
663 while ((chunk = (THnSparseArrayChunk*) iChunk())) {
664 const Int_t chunkSize = chunk->GetEntries();
665 Char_t* buf = chunk->fCoordinates;
666 const Int_t singleCoordSize = chunk->fSingleCoordinateSize;
667 const Char_t* endbuf = buf + singleCoordSize * chunkSize;
668 for (; buf < endbuf; buf += singleCoordSize, ++idx) {
669 Long64_t hash = compactCoord.GetHashFromBuffer(buf);
670 Long64_t linidx = fBins.GetValue(hash);
671 if (linidx) {
672 Long64_t nextidx = linidx;
673 while (nextidx) {
674 // must be a collision, so go to fBinsContinued.
675 linidx = nextidx;
676 nextidx = fBinsContinued.GetValue(linidx);
677 }
678 fBinsContinued.Add(linidx, idx + 1);
679 } else {
680 fBins.Add(hash, idx + 1);
681 }
682 }
683 }
684}
685
686////////////////////////////////////////////////////////////////////////////////
687/// Initialize storage for nbins
688
690 if (!fBins.GetSize() && fBinContent.GetSize()) {
691 FillExMap();
692 }
693 if (2 * nbins > fBins.Capacity()) {
694 fBins.Expand(3 * nbins);
695 }
696}
697
698////////////////////////////////////////////////////////////////////////////////
699/// Get the bin index for the n dimensional tuple x,
700/// allocate one if it doesn't exist yet and "allocate" is true.
701
702Long64_t THnSparse::GetBin(const Double_t* x, Bool_t allocate /* = kTRUE */)
703{
704 THnSparseCompactBinCoord* cc = GetCompactCoord();
705 Int_t *coord = cc->GetCoord();
706 for (Int_t i = 0; i < fNdimensions; ++i)
707 coord[i] = GetAxis(i)->FindBin(x[i]);
708 cc->UpdateCoord();
709
710 return GetBinIndexForCurrentBin(allocate);
711}
712
713
714////////////////////////////////////////////////////////////////////////////////
715/// Get the bin index for the n dimensional tuple addressed by "name",
716/// allocate one if it doesn't exist yet and "allocate" is true.
717
718Long64_t THnSparse::GetBin(const char* name[], Bool_t allocate /* = kTRUE */)
719{
720 THnSparseCompactBinCoord* cc = GetCompactCoord();
721 Int_t *coord = cc->GetCoord();
722 for (Int_t i = 0; i < fNdimensions; ++i)
723 coord[i] = GetAxis(i)->FindBin(name[i]);
724 cc->UpdateCoord();
725
726 return GetBinIndexForCurrentBin(allocate);
727}
728
729////////////////////////////////////////////////////////////////////////////////
730/// Get the bin index for the n dimensional coordinates coord,
731/// allocate one if it doesn't exist yet and "allocate" is true.
732
733Long64_t THnSparse::GetBin(const Int_t* coord, Bool_t allocate /*= kTRUE*/)
734{
735 GetCompactCoord()->SetCoord(coord);
736 return GetBinIndexForCurrentBin(allocate);
737}
738
739////////////////////////////////////////////////////////////////////////////////
740/// Return the content of the filled bin number "idx".
741/// If coord is non-null, it will contain the bin's coordinates for each axis
742/// that correspond to the bin.
743
745{
746 if (idx >= 0) {
748 idx %= fChunkSize;
749 if (chunk && chunk->fContent->GetSize() > idx) {
750 if (coord) {
751 THnSparseCompactBinCoord* cc = GetCompactCoord();
752 Int_t sizeCompact = cc->GetBufferSize();
753 cc->SetCoordFromBuffer(chunk->fCoordinates + idx * sizeCompact,
754 coord);
755
756 }
757 return chunk->fContent->GetAt(idx);
758 }
759 }
760 if (coord)
761 memset(coord, -1, sizeof(Int_t) * fNdimensions);
762 return 0.;
763}
764
765////////////////////////////////////////////////////////////////////////////////
766/// Get square of the error of bin addressed by linidx as
767/// \f$\sum weight^{2}\f$
768/// If errors are not enabled (via Sumw2() or CalculateErrors())
769/// return contents.
770
772 if (!GetCalculateErrors())
773 return GetBinContent(linidx);
774
775 if (linidx < 0) return 0.;
776 THnSparseArrayChunk* chunk = GetChunk(linidx / fChunkSize);
777 linidx %= fChunkSize;
778 if (!chunk || chunk->fContent->GetSize() < linidx)
779 return 0.;
780
781 return chunk->fSumw2->GetAt(linidx);
782}
783
784
785////////////////////////////////////////////////////////////////////////////////
786/// Return the index for fCurrentBinIndex.
787/// If it doesn't exist then return -1, or allocate a new bin if allocate is set
788
790{
791 THnSparseCompactBinCoord* cc = GetCompactCoord();
792 ULong64_t hash = cc->GetHash();
793 if (fBinContent.GetSize() && !fBins.GetSize())
794 FillExMap();
795 Long64_t linidx = (Long64_t) fBins.GetValue(hash);
796 while (linidx) {
797 // fBins stores index + 1!
798 THnSparseArrayChunk* chunk = GetChunk((linidx - 1)/ fChunkSize);
799 if (chunk->Matches((linidx - 1) % fChunkSize, cc->GetBuffer()))
800 return linidx - 1; // we store idx+1, 0 is "TExMap: not found"
801
802 Long64_t nextlinidx = fBinsContinued.GetValue(linidx);
803 if (!nextlinidx) break;
804
805 linidx = nextlinidx;
806 }
807 if (!allocate) return -1;
808
809 ++fFilledBins;
810
811 // allocate bin in chunk
813 Long64_t newidx = chunk ? ((Long64_t) chunk->GetEntries()) : -1;
814 if (!chunk || newidx == (Long64_t)fChunkSize) {
815 chunk = AddChunk();
816 newidx = 0;
817 }
818 chunk->AddBin(newidx, cc->GetBuffer());
819
820 // store translation between hash and bin
821 newidx += (fBinContent.GetEntriesFast() - 1) * fChunkSize;
822 if (!linidx) {
823 // fBins didn't find it
824 if (2 * GetNbins() > fBins.Capacity())
825 fBins.Expand(3 * GetNbins());
826 fBins.Add(hash, newidx + 1);
827 } else {
828 // fBins contains one, but it's the wrong one;
829 // add entry to fBinsContinued.
830 fBinsContinued.Add(linidx, newidx + 1);
831 }
832 return newidx;
833}
834
835////////////////////////////////////////////////////////////////////////////////
836/// Return THnSparseCompactBinCoord object.
837
838THnSparseCompactBinCoord* THnSparse::GetCompactCoord() const
839{
840 if (!fCompactCoord) {
841 Int_t *bins = new Int_t[fNdimensions];
842 for (Int_t d = 0; d < fNdimensions; ++d)
843 bins[d] = GetAxis(d)->GetNbins();
844 const_cast<THnSparse*>(this)->fCompactCoord
845 = new THnSparseCompactBinCoord(fNdimensions, bins);
846 delete [] bins;
847 }
848 return fCompactCoord;
849}
850
851////////////////////////////////////////////////////////////////////////////////
852/// Return the amount of filled bins over all bins
853
855 Double_t nbinsTotal = 1.;
856 for (Int_t d = 0; d < fNdimensions; ++d)
857 nbinsTotal *= GetAxis(d)->GetNbins() + 2;
858 return fFilledBins / nbinsTotal;
859}
860
861////////////////////////////////////////////////////////////////////////////////
862/// Return the amount of used memory over memory that would be used by a
863/// non-sparse n-dimensional histogram. The value is approximate.
864
866 Int_t arrayElementSize = 0;
867 if (fFilledBins) {
868 TClass* clArray = GetChunk(0)->fContent->IsA();
869 TDataMember* dm = clArray ? clArray->GetDataMember("fArray") : 0;
870 arrayElementSize = dm ? dm->GetDataType()->Size() : 0;
871 }
872 if (!arrayElementSize) {
873 Warning("GetSparseFractionMem", "Cannot determine type of elements!");
874 return -1.;
875 }
876
877 Double_t sizePerChunkElement = arrayElementSize + GetCompactCoord()->GetBufferSize();
878 if (fFilledBins && GetChunk(0)->fSumw2)
879 sizePerChunkElement += sizeof(Double_t); /* fSumw2 */
880
881 Double_t size = 0.;
882 size += fBinContent.GetEntries() * (GetChunkSize() * sizePerChunkElement + sizeof(THnSparseArrayChunk));
883 size += + 3 * sizeof(Long64_t) * fBins.GetSize() /* TExMap */;
884
885 Double_t nbinsTotal = 1.;
886 for (Int_t d = 0; d < fNdimensions; ++d)
887 nbinsTotal *= GetAxis(d)->GetNbins() + 2;
888
889 return size / nbinsTotal / arrayElementSize;
890}
891
892////////////////////////////////////////////////////////////////////////////////
893/// Create an iterator over all filled bins of a THnSparse.
894/// Use THnIter instead.
895
897{
898 return new THnSparseBinIter(respectAxisRange, this);
899}
900
901////////////////////////////////////////////////////////////////////////////////
902/// Set content of bin with index "bin" to "v"
903
905{
907 chunk->fContent->SetAt(v, bin % fChunkSize);
908 ++fEntries;
909}
910
911////////////////////////////////////////////////////////////////////////////////
912/// Set error of bin with index "bin" to "e", enable errors if needed
913
915{
917 if (!chunk->fSumw2 ) {
918 // if fSumw2 is zero GetCalculateErrors should return false
919 if (GetCalculateErrors()) {
920 Error("SetBinError", "GetCalculateErrors() logic error!");
921 }
922 Sumw2(); // enable error calculation
923 }
924
925 chunk->fSumw2->SetAt(e2, bin % fChunkSize);
926}
927
928////////////////////////////////////////////////////////////////////////////////
929/// Add "e" to error of bin with index "bin", enable errors if needed
930
932{
934 if (!chunk->fSumw2 ) {
935 // if fSumw2 is zero GetCalculateErrors should return false
936 if (GetCalculateErrors()) {
937 Error("SetBinError", "GetCalculateErrors() logic error!");
938 }
939 Sumw2(); // enable error calculation
940 }
941
942 (*chunk->fSumw2)[bin % fChunkSize] += e2;
943}
944
945////////////////////////////////////////////////////////////////////////////////
946/// Enable calculation of errors
947
949{
950 if (GetCalculateErrors()) return;
951
952 fTsumw2 = 0.;
953 TIter iChunk(&fBinContent);
954 THnSparseArrayChunk* chunk = 0;
955 while ((chunk = (THnSparseArrayChunk*) iChunk()))
956 chunk->Sumw2();
957}
958
959////////////////////////////////////////////////////////////////////////////////
960/// Clear the histogram
961
962void THnSparse::Reset(Option_t *option /*= ""*/)
963{
964 fFilledBins = 0;
965 fBins.Delete();
968 ResetBase(option);
969}
970
ROOT::R::TRInterface & r
Definition Object.C:4
#define d(i)
Definition RSha256.hxx:102
int Int_t
Definition RtypesCore.h:45
unsigned char UChar_t
Definition RtypesCore.h:38
char Char_t
Definition RtypesCore.h:37
unsigned int UInt_t
Definition RtypesCore.h:46
bool Bool_t
Definition RtypesCore.h:63
double Double_t
Definition RtypesCore.h:59
long long Long64_t
Definition RtypesCore.h:73
unsigned long long ULong64_t
Definition RtypesCore.h:74
const Bool_t kTRUE
Definition RtypesCore.h:91
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:364
char name[80]
Definition TGX11.cxx:110
float xmin
float xmax
Binding & operator=(OUT(*fun)(void))
Iterator over THnBase bins (internal implementation).
Definition THnBase.h:285
virtual Int_t GetCoord(Int_t dim) const =0
virtual Long64_t Next(Int_t *coord=0)=0
Array of doubles (64 bits per element).
Definition TArrayD.h:27
void SetAt(Double_t v, Int_t i)
Definition TArrayD.h:51
Double_t GetAt(Int_t i) const
Definition TArrayD.h:45
Double_t * fArray
Definition TArrayD.h:30
Abstract array base class.
Definition TArray.h:31
virtual void SetAt(Double_t v, Int_t i)=0
virtual Double_t GetAt(Int_t i) const =0
Int_t GetSize() const
Definition TArray.h:47
virtual Int_t FindBin(Double_t x)
Find bin number corresponding to abscissa x.
Definition TAxis.cxx:293
Int_t GetNbins() const
Definition TAxis.h:121
TClass instances represent classes, structs and namespaces in the ROOT type system.
Definition TClass.h:80
TDataMember * GetDataMember(const char *datamember) const
Return pointer to datamember object with name "datamember".
Definition TClass.cxx:3416
virtual void SetOwner(Bool_t enable=kTRUE)
Set whether this collection is the owner (enable==true) of its content.
virtual Int_t GetSize() const
Return the capacity of the collection, i.e.
All ROOT classes may have RTTI (run time type identification) support added.
Definition TDataMember.h:31
TDataType * GetDataType() const
Definition TDataMember.h:76
Int_t Size() const
Get size of basic typedef'ed type.
void Delete(Option_t *opt="")
Delete all entries stored in the TExMap.
Definition TExMap.cxx:163
void Expand(Int_t newsize)
Expand the TExMap.
Definition TExMap.cxx:278
Int_t GetSize() const
Definition TExMap.h:71
void Add(ULong64_t hash, Long64_t key, Long64_t value)
Add an (key,value) pair to the table. The key should be unique.
Definition TExMap.cxx:87
Long64_t GetValue(ULong64_t hash, Long64_t key)
Return the value belonging to specified key and hash value.
Definition TExMap.cxx:173
Int_t Capacity() const
Definition TExMap.h:69
Multidimensional histogram base.
Definition THnBase.h:43
Double_t fEntries
browser-helpers for each axis
Definition THnBase.h:48
Int_t GetNdimensions() const
Definition THnBase.h:135
void ResetBase(Option_t *option="")
Clear the histogram.
Definition THnBase.cxx:1151
Bool_t GetCalculateErrors() const
Definition THnBase.h:136
Double_t fTsumw2
Definition THnBase.h:50
TAxis * GetAxis(Int_t dim) const
Definition THnBase.h:125
Int_t fNdimensions
Definition THnBase.h:45
THnSparseArrayChunk is used internally by THnSparse.
Bool_t Matches(Int_t idx, const Char_t *idxbuf) const
virtual ~THnSparseArrayChunk()
Destructor.
Int_t fSingleCoordinateSize
size of the allocated coordinate buffer; -1 means none or fCoordinatesSize
void Sumw2()
Turn on support of errors.
void AddBin(Int_t idx, const Char_t *idxbuf)
Create a new bin in this chunk.
Efficient multidimensional histogram.
Definition THnSparse.h:36
Double_t GetSparseFractionBins() const
Return the amount of filled bins over all bins.
Long64_t GetBin(const Int_t *idx) const
Definition THnSparse.h:95
Double_t GetSparseFractionMem() const
Return the amount of used memory over memory that would be used by a non-sparse n-dimensional histogr...
void SetBinContent(const Int_t *idx, Double_t v)
Forwards to THnBase::SetBinContent().
Definition THnSparse.h:104
void InitStorage(Int_t *nbins, Int_t chunkSize)
Initialize the storage of a histogram created via Init()
THnSparseArrayChunk * GetChunk(Int_t idx) const
Definition THnSparse.h:55
THnSparseCompactBinCoord * fCompactCoord
filled bins for non-unique hashes, containing pairs of (bin index 0, bin index 1)
Definition THnSparse.h:43
Int_t GetChunkSize() const
Definition THnSparse.h:87
virtual ~THnSparse()
Destruct a THnSparse.
ROOT::Internal::THnBaseBinIter * CreateIter(Bool_t respectAxisRange) const
Create an iterator over all filled bins of a THnSparse.
Double_t GetBinContent(const Int_t *idx) const
Forwards to THnBase::GetBinContent() overload.
Definition THnSparse.h:120
Long64_t fFilledBins
Definition THnSparse.h:39
void Reset(Option_t *option="")
Clear the histogram.
TObjArray fBinContent
Definition THnSparse.h:40
THnSparseCompactBinCoord * GetCompactCoord() const
Return THnSparseCompactBinCoord object.
Double_t GetBinError2(Long64_t linidx) const
Get square of the error of bin addressed by linidx as If errors are not enabled (via Sumw2() or Calc...
THnSparse()
Construct an empty THnSparse.
THnSparseArrayChunk * AddChunk()
Create a new chunk of bin content.
virtual TArray * GenerateArray() const =0
void FillExMap()
We have been streamed; set up fBins.
void AddBinError2(Long64_t bin, Double_t e2)
Add "e" to error of bin with index "bin", enable errors if needed.
TExMap fBins
Definition THnSparse.h:41
TExMap fBinsContinued
filled bins
Definition THnSparse.h:42
void Sumw2()
Enable calculation of errors.
void Reserve(Long64_t nbins)
Initialize storage for nbins.
void AddBinContent(const Int_t *idx, Double_t v=1.)
Forwards to THnBase::SetBinContent().
Definition THnSparse.h:112
void SetBinError2(Long64_t bin, Double_t e2)
Set error of bin with index "bin" to "e", enable errors if needed.
Int_t fChunkSize
Definition THnSparse.h:38
Long64_t GetNbins() const
Definition THnSparse.h:92
Long64_t GetBinIndexForCurrentBin(Bool_t allocate)
Return the index for fCurrentBinIndex.
Int_t GetEntriesFast() const
Definition TObjArray.h:64
virtual void AddLast(TObject *obj)
Add object in the next empty slot in the array.
TObject * Last() const
Return the object in the last filled slot. Returns 0 if no entries.
Int_t GetEntries() const
Return the number of objects in array (i.e.
virtual void Delete(Option_t *option="")
Remove all objects from the array AND delete all heap based objects.
virtual void Clear(Option_t *="")
Definition TObject.h:115
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:879
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:893
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
Py_ssize_t GetBuffer(PyObject *pyobject, char tc, int size, void *&buf, bool check=true)
Definition Utility.cxx:614
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...