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