Logo ROOT   6.12/07
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 
20 namespace {
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 
51 Int_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 
59 Long64_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
100 THnSparseCoordCompression is a class used by THnSparse internally. It
101 represents a compacted n-dimensional array of bin coordinates (indices).
102 As the total number of bins in each dimension is known by THnSparse, bin
103 indices can be compacted to only use the amount of bins needed by the total
104 number 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
106 28 bits (4+7+1+5+4+7), i.e. less than a 32bit integer. The tricky part is
107 the fast compression and decompression, the platform-independent storage
108 (think of endianness: the bits of the number 0x123456 depend on the
109 platform), and the hashing needed by THnSparseArrayChunk.
110 */
111 
112 
113 class THnSparseCoordCompression {
114 public:
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 
127 protected:
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  }
134 private:
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 
150 THnSparseCoordCompression::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 
168 THnSparseCoordCompression::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 
180 THnSparseCoordCompression& 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 
196 THnSparseCoordCompression::~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 
206 void 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 
234 ULong64_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 
270 ULong64_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 
318 ULong64_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
347 THnSparseCompactBinCoord is a class used by THnSparse internally. It
348 maps between an n-dimensional array of bin coordinates (indices) and
349 its compact version, the THnSparseCoordCompression.
350 */
351 
352 class THnSparseCompactBinCoord: public THnSparseCoordCompression {
353 public:
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 
371 private:
372  // intentionally not implemented
373  THnSparseCompactBinCoord(const THnSparseCompactBinCoord&);
374  // intentionally not implemented
375  THnSparseCompactBinCoord& operator=(const THnSparseCompactBinCoord&);
376 
377 private:
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 
392 THnSparseCompactBinCoord::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 
407 THnSparseCompactBinCoord::~THnSparseCompactBinCoord()
408 {
409  delete [] fCoordBuffer;
410  delete [] fCurrentBin;
411 }
412 
413 /** \class THnSparseArrayChunk
414 THnSparseArrayChunk is used internally by THnSparse.
415 THnSparse stores its (dynamic size) array of bin coordinates and their
416 contents (and possibly errors) in a TObjArray of THnSparseArrayChunk. Each
417 of the chunks holds an array of THnSparseCompactBinCoord and the content
418 (a TArray*), which is created outside (by the templated derived classes of
419 THnSparse) 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 
451 void 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) {
463  Int_t chunksize = fSingleCoordinateSize * fContent->GetSize();
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 
474  memcpy(fCoordinates + idx * fSingleCoordinateSize, coordbuf, fSingleCoordinateSize);
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 
496 Efficient multidimensional histogram.
497 
498 Use a THnSparse instead of TH1 / TH2 / TH3 / array for histogramming when
499 only a small fraction of bins is filled. A 10-dimensional histogram with 10
500 bins per dimension has 10^10 bins; in a naive implementation this will not
501 fit in memory. THnSparse only allocates memory for the bins that have
502 non-zero bin content instead, drastically reducing both the memory usage
503 and the access time.
504 
505 To construct a THnSparse object you must use one of its templated, derived
506 classes:
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 
514 They take name and title, the number of dimensions, and for each dimension
515 the number of bins, the minimal, and the maximal value on the dimension's
516 axis. 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
524 A THnSparse is filled just like a regular histogram, using
525 THnSparse::Fill(x, weight), where x is a n-dimensional Double_t value.
526 To take errors into account, Sumw2() must be called before filling the
527 histogram.
528 
529 Bins are allocated as needed; the status of the allocation can be observed
530 by GetSparseFractionBins(), GetSparseFractionMem().
531 
532 ## Fast Bin Content Access
533 When iterating over a THnSparse one should only look at filled bins to save
534 processing time. The number of filled bins is returned by
535 THnSparse::GetNbins(); the bin content for each (linear) bin number can
536 be retrieved by THnSparse::GetBinContent(linidx, (Int_t*)coord).
537 After the call, coord will contain the bin coordinate of each axis for the bin
538 with 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
545 TH1 and TH2 are generally faster than THnSparse for one and two dimensional
546 distributions. THnSparse becomes competitive for a sparsely filled TH3
547 with large numbers of bins per dimension. The tutorial sparsehist.C
548 shows the turning point. On a AMD64 with 8GB memory, THnSparse "wins"
549 starting with a TH3 with 30 bins per dimension. Using a THnSparse for a
550 one-dimensional histogram is only reasonable if it has a huge number of bins.
551 
552 ## Projections
553 The dimensionality of a THnSparse can be reduced by projecting it to
554 1, 2, 3, or n dimensions, which can be represented by a TH1, TH2, TH3, or
555 a THnSparse. See the Projection() members. To only project parts of the
556 histogram, call
557 
558  THnSparse::GetAxis(12)->SetRange(from_bin, to_bin);
559 
560 ## Internal Representation
561 An entry for a filled bin consists of its n-dimensional coordinates and
562 its bin content. The coordinates are compacted to use as few bits as
563 possible; e.g. a histogram with 10 bins in x and 20 bins in y will only
564 use 4 bits for the x representation and 5 bits for the y representation.
565 This is handled by the internal class THnSparseCompactBinCoord.
566 Bin data (content and coordinates) are allocated in chunks of size
567 fChunkSize; this parameter can be set when constructing a THnSparse. Each
568 chunk is represented by an object of class THnSparseArrayChunk.
569 
570 Translation from an n-dimensional bin coordinate to the linear index within
571 the chunks is done by GetBin(). It creates a hash from the compacted bin
572 coordinates (the hash of a bin coordinate is the compacted coordinate itself
573 if it takes less than 8 bytes, the size of a Long64_t.
574 This hash is used to lookup the linear index in the TExMap member fBins;
575 the coordinates of the entry fBins points to is compared to the coordinates
576 passed to GetBin(). If they do not match, these two coordinates have the same
577 hash - which is extremely unlikely but (for the case where the compact bin
578 coordinates are larger than 4 bytes) possible. In this case, fBinsContinued
579 contains a chain of linear indexes with the same hash. Iterating through this
580 chain and comparing each bin coordinates with the one passed to GetBin() will
581 retrieve 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 
604 THnSparse::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 {
626  THnSparseArrayChunk* chunk = GetChunk(bin / fChunkSize);
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 
647 void 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 
703 Long64_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 
719 Long64_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 
734 Long64_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 
745 Double_t THnSparse::GetBinContent(Long64_t idx, Int_t* coord /* = 0 */) const
746 {
747  if (idx >= 0) {
748  THnSparseArrayChunk* chunk = GetChunk(idx / fChunkSize);
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 
839 THnSparseCompactBinCoord* 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 {
907  THnSparseArrayChunk* chunk = GetChunk(bin / fChunkSize);
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 {
917  THnSparseArrayChunk* chunk = GetChunk(bin / fChunkSize);
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 {
934  THnSparseArrayChunk* chunk = GetChunk(bin / fChunkSize);
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 
963 void THnSparse::Reset(Option_t *option /*= ""*/)
964 {
965  fFilledBins = 0;
966  fBins.Delete();
969  ResetBase(option);
970 }
971 
Abstract array base class.
Definition: TArray.h:31
virtual void Clear(Option_t *="")
Definition: TObject.h:100
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
THnSparseArrayChunk is used internally by THnSparse.
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
void AddBinContent(const Int_t *idx, Double_t v=1.)
Definition: THnSparse.h:108
float xmin
Definition: THbookFile.cxx:93
Bool_t GetCalculateErrors() const
Definition: THnBase.h:136
void SetAt(Double_t v, Int_t i)
Definition: TArrayD.h:51
long long Long64_t
Definition: RtypesCore.h:69
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
void FillExMap()
We have been streamed; set up fBins.
Definition: THnSparse.cxx:656
Double_t fTsumw2
Definition: THnBase.h:50
Int_t GetSize() const
Definition: TExMap.h:71
const char Option_t
Definition: RtypesCore.h:62
THnSparseArrayChunk * AddChunk()
Create a new chunk of bin content.
Definition: THnSparse.cxx:635
virtual void Delete(Option_t *option="")
Remove all objects from the array AND delete all heap based objects.
Definition: TObjArray.cxx:355
All ROOT classes may have RTTI (run time type identification) support added.
Definition: TDataMember.h:31
virtual ~THnSparseArrayChunk()
Destructor.
Definition: THnSparse.cxx:441
Long64_t GetBin(const Int_t *idx) const
Definition: THnSparse.h:94
virtual void SetOwner(Bool_t enable=kTRUE)
Set whether this collection is the owner (enable==true) of its content.
Long64_t fFilledBins
Definition: THnSparse.h:39
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
void Sumw2()
Turn on support of errors.
Definition: THnSparse.cxx:481
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
Int_t GetEntries() const
void InitStorage(Int_t *nbins, Int_t chunkSize)
Initialize the storage of a histogram created via Init()
Definition: THnSparse.cxx:647
THnSparseCompactBinCoord * fCompactCoord
filled bins for non-unique hashes, containing pairs of (bin index 0, bin index 1) ...
Definition: THnSparse.h:43
Long64_t GetNbins() const
Definition: THnSparse.h:91
Double_t fEntries
browser-helpers for each axis
Definition: THnBase.h:48
TObject * Last() const
Return the object in the last filled slot. Returns 0 if no entries.
Definition: TObjArray.cxx:505
Bool_t Matches(Int_t idx, const Char_t *idxbuf) const
Double_t x[n]
Definition: legend1.C:17
virtual void SetAt(Double_t v, Int_t i)=0
Efficient multidimensional histogram.
Definition: THnSparse.h:36
virtual TArray * GenerateArray() const =0
Double_t * fArray
Definition: TArrayD.h:30
TDataType * GetDataType() const
Definition: TDataMember.h:74
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition: TObject.h:271
Int_t fSingleCoordinateSize
size of the allocated coordinate buffer; -1 means none or fCoordinatesSize
THnSparseArrayChunk * GetChunk(Int_t idx) const
Definition: THnSparse.h:55
virtual Int_t GetCoord(Int_t dim) const =0
virtual ~THnSparse()
Destruct a THnSparse.
Definition: THnSparse.cxx:617
Int_t Capacity() const
Definition: TExMap.h:69
Double_t GetBinContent(const Int_t *idx) const
Definition: THnSparse.h:116
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
TExMap fBinsContinued
filled bins
Definition: THnSparse.h:42
void Reset(Option_t *option="")
Clear the histogram.
Definition: THnSparse.cxx:963
TAxis * GetAxis(Int_t dim) const
Definition: THnBase.h:125
Long64_t GetValue(ULong64_t hash, Long64_t key)
Return the value belonging to specified key and hash value.
Definition: TExMap.cxx:173
ROOT::R::TRInterface & r
Definition: Object.C:4
TExMap fBins
Definition: THnSparse.h:41
Int_t GetSize() const
Definition: TArray.h:47
SVector< double, 2 > v
Definition: Dict.h:5
Iterator over THnBase bins (internal implementation).
Definition: THnBase.h:286
void Sumw2()
Enable calculation of errors.
Definition: THnSparse.cxx:949
unsigned int UInt_t
Definition: RtypesCore.h:42
Int_t GetEntriesFast() const
Definition: TObjArray.h:64
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
The ROOT global object gROOT contains a list of all defined classes.
Definition: TClass.h:75
float xmax
Definition: THbookFile.cxx:93
void SetBinContent(const Int_t *idx, Double_t v)
Definition: THnSparse.h:101
void Delete(Option_t *opt="")
Delete all entries stored in the TExMap.
Definition: TExMap.cxx:163
ROOT::Internal::THnBaseBinIter * CreateIter(Bool_t respectAxisRange) const
Create an iterator over all filled bins of a THnSparse.
Definition: THnSparse.cxx:897
virtual Int_t FindBin(Double_t x)
Find bin number corresponding to abscissa x.
Definition: TAxis.cxx:279
void Reserve(Long64_t nbins)
Initialize storage for nbins.
Definition: THnSparse.cxx:690
#define ClassImp(name)
Definition: Rtypes.h:359
double Double_t
Definition: RtypesCore.h:55
unsigned long long ULong64_t
Definition: RtypesCore.h:70
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
Double_t GetSparseFractionBins() const
Return the amount of filled bins over all bins.
Definition: THnSparse.cxx:855
THnSparse()
Construct an empty THnSparse.
Definition: THnSparse.cxx:590
Array of doubles (64 bits per element).
Definition: TArrayD.h:27
Int_t GetNdimensions() const
Definition: THnBase.h:135
Binding & operator=(OUT(*fun)(void))
void AddBin(Int_t idx, const Char_t *idxbuf)
Create a new bin in this chunk.
Definition: THnSparse.cxx:451
char Char_t
Definition: RtypesCore.h:29
THnSparseCompactBinCoord * GetCompactCoord() const
Return THnSparseCompactBinCoord object.
Definition: THnSparse.cxx:839
virtual Double_t GetAt(Int_t i) const =0
Int_t fNdimensions
Definition: THnBase.h:45
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:539
Int_t Size() const
Get size of basic typedef&#39;ed type.
Definition: TDataType.cxx:366
Bool_t RespectsAxisRange() const
Definition: THnBase.h:292
Double_t GetAt(Int_t i) const
Definition: TArrayD.h:45
Int_t fChunkSize
Definition: THnSparse.h:38
Int_t GetEntries() const
Return the number of objects in array (i.e.
Definition: TObjArray.cxx:522
virtual void AddLast(TObject *obj)
Add object in the next empty slot in the array.
Definition: TObjArray.cxx:177
Multidimensional histogram base.
Definition: THnBase.h:43
TDataMember * GetDataMember(const char *datamember) const
Return pointer to datamember object with name "datamember".
Definition: TClass.cxx:3280
unsigned char UChar_t
Definition: RtypesCore.h:34
Int_t GetNbins() const
Definition: TAxis.h:121
virtual Int_t GetSize() const
Definition: TCollection.h:180
Int_t GetChunkSize() const
Definition: THnSparse.h:86
const Bool_t kTRUE
Definition: RtypesCore.h:87
Long64_t GetBinIndexForCurrentBin(Bool_t allocate)
Return the index for fCurrentBinIndex.
Definition: THnSparse.cxx:790
const Int_t n
Definition: legend1.C:16
virtual Long64_t Next(Int_t *coord=0)=0
TObjArray fBinContent
Definition: THnSparse.h:40
char name[80]
Definition: TGX11.cxx:109
void ResetBase(Option_t *option="")
Clear the histogram.
Definition: THnBase.cxx:1151
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:866
void Expand(Int_t newsize)
Expand the TExMap.
Definition: TExMap.cxx:278