reduce size on disk for storing PDF nodes
authorabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 16 Dec 2009 12:02:43 +0000 (12:02 +0000)
committerabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 16 Dec 2009 12:02:43 +0000 (12:02 +0000)
 - go from TVectoT to Double_t[]
 - store selectively interpolation covariance matrix
 - bootstrap data members which are not saved on disk

STAT/TKDInterpolator.cxx
STAT/TKDInterpolator.h
STAT/TKDInterpolatorBase.cxx
STAT/TKDInterpolatorBase.h
STAT/TKDNodeInfo.cxx
STAT/TKDNodeInfo.h
STAT/TKDPDF.cxx
STAT/TKDPDF.h
STAT/TKDSpline.cxx
STAT/TKDSpline.h

index 8c4ddd5..31d846a 100644 (file)
@@ -34,19 +34,20 @@ TKDInterpolator::~TKDInterpolator()
 //_________________________________________________________________
 void TKDInterpolator::AddNode(const TKDNodeInfo &node)
 {
-  if(!fTNodes){
+  if(!fNodes){
     Warning("TKDInterpolator::SetNode()", "Node array not defined.");
     return;
   }
 
-  new((*fTNodes)[fNTNodes++]) TKDNodeInfo(node);
+  Int_t n(GetNTNodes());
+  new((*fNodes)[n++]) TKDNodeInfo(node);
 }
 
 //_________________________________________________________________
-void TKDInterpolator::Build(Int_t npoints, Int_t ndim)
+Bool_t TKDInterpolator::Build(Int_t npoints, Int_t ndim)
 {
   fNSize = ndim;
-  TKDInterpolatorBase::Build(npoints);
+  return TKDInterpolatorBase::Build(npoints);
 }
 
 //_________________________________________________________________
@@ -57,8 +58,8 @@ Int_t TKDInterpolator::GetNodeIndex(const Float_t *p)
 //   for(int i=0; i<fNSize; i++) printf("%f ", p[i]);
 //   printf("] ...\n");
 
-  for(Int_t i=fNTNodes; i--;)
-    if(((TKDNodeInfo*)(*fTNodes)[i])->Has(p)) return i;
+  for(Int_t i=GetNTNodes(); i--;)
+    if(((TKDNodeInfo*)(*fNodes)[i])->Has(p)) return i;
   
   printf("Point p[");
   for(int i=0; i<fNSize; i++) printf("%f ", p[i]);
@@ -70,15 +71,15 @@ Int_t TKDInterpolator::GetNodeIndex(const Float_t *p)
 //_________________________________________________________________
 Bool_t TKDInterpolator::SetNode(Int_t inode, const TKDNodeInfo &ref)
 {
-  if(!fTNodes){
+  if(!fNodes){
     Warning("TKDInterpolator::SetNode()", "Node array not defined.");
     return kFALSE;
   }
-  if(inode >= fNTNodes){
-    Warning("TKDInterpolator::SetNode()", Form("Node array defined up to %d.", fNTNodes));
+  if(inode >= GetNTNodes()){
+    Warning("TKDInterpolator::SetNode()", Form("Node array defined up to %d.", GetNTNodes()));
     return kFALSE;
   }
-  TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode];
+  TKDNodeInfo *node = (TKDNodeInfo*)(*fNodes)[inode];
   (*node) = ref;
   return kTRUE;
 }
index 1c3f6a9..73b09e4 100644 (file)
@@ -12,8 +12,8 @@ public:
   TKDInterpolator(Int_t ndim, Int_t npoints=0);
   ~TKDInterpolator();
   void       AddNode(const TKDNodeInfo &ref);
-  void       Build(Int_t ndim) {TKDInterpolatorBase::Build(ndim);}
-  void       Build(Int_t npoints, Int_t ndim);
+  Bool_t     Build(Int_t ndim) {return TKDInterpolatorBase::Build(ndim);}
+  Bool_t     Build(Int_t npoints, Int_t ndim);
   Int_t      GetNodeIndex(const Float_t *p);
   Bool_t     SetNode(Int_t i, const TKDNodeInfo &ref);
 
index 074c434..ff847c9 100644 (file)
@@ -30,17 +30,16 @@ ClassImp(TKDInterpolatorBase)
 //_________________________________________________________________
 TKDInterpolatorBase::TKDInterpolatorBase(Int_t dim) :
   fNSize(dim)
-  ,fNTNodes(0)
-  ,fTNodes(0x0)
-  ,fTNodesDraw(0x0)
+  ,fNodes(NULL)
+  ,fNodesDraw(NULL)
   ,fStatus(0)
   ,fLambda(1 + dim + (dim*(dim+1)>>1))
   ,fDepth(-1)
   ,fAlpha(.5)
-  ,fRefPoints(0x0)
-  ,fBuffer(0x0)
-  ,fKDhelper(0x0)
-  ,fFitter(0x0)
+  ,fRefPoints(NULL)
+  ,fBuffer(NULL)
+  ,fKDhelper(NULL)
+  ,fFitter(NULL)
 {
 // Default constructor. To be used with care since in this case building
 // of data structure is completly left to the user responsability.
@@ -48,18 +47,46 @@ TKDInterpolatorBase::TKDInterpolatorBase(Int_t dim) :
 }
 
 //_________________________________________________________________
-void   TKDInterpolatorBase::Build(Int_t n)
+Bool_t TKDInterpolatorBase::Build(Int_t n)
 {
   // allocate memory for data
+  if(Int_t((1.+fAlpha)*fLambda) > n){ // check granularity
+    Error("TKDInterpolatorBase::Build()", Form("Minimum number of points [%d] needed for interpolation exceeds number of evaluation points [%d]. Please increase granularity.", Int_t((1.+fAlpha)*fLambda), n));
+    return kFALSE;
+  }
 
-  if(fTNodes) delete fTNodes;
-  fNTNodes = n;
-  // check granularity
-  if(Int_t((1.+fAlpha)*fLambda) > fNTNodes){
-    Warning("TKDInterpolatorBase::Build()", Form("Minimum number of points [%d] needed for interpolation exceeds number of evaluation points [%d]. Please increase granularity.", Int_t((1.+fAlpha)*fLambda), fNTNodes));
+  if(fNodes){
+    Warning("TKDInterpolatorBase::Build()", "Data already allocated.");
+    fNodes->Delete();
+  } else {
+    fNodes = new TClonesArray("TKDNodeInfo", n); 
+    fNodes->SetOwner();
   }
-  fTNodes = new TClonesArray("TKDNodeInfo", fNTNodes);
-  for(int in=0; in<fNTNodes; in++) new ((*fTNodes)[in]) TKDNodeInfo(fNSize);
+
+  for(int in=0; in<n; in++) new ((*fNodes)[in]) TKDNodeInfo(fNSize);
+
+  return kTRUE;
+}
+
+//_________________________________________________________________
+Bool_t TKDInterpolatorBase::Bootstrap()
+{
+  if(!fNodes){
+    Error("TKDInterpolatorBase::Bootstrap()", "Nodes missing. Nothing to bootstrap from.");
+    return kFALSE;
+  }
+  Int_t in = GetNTNodes(); TKDNodeInfo *n(NULL);
+  while(in--){ 
+    if(!(n=(TKDNodeInfo*)(*fNodes)[in])){
+      Error("TKDInterpolatorBase::Bootstrap()", Form("Node @ %d missing.", in));
+      return kFALSE;
+    }
+    n->Bootstrap();
+    if(!fNSize) fNSize  = n->GetDimension();
+    //n->SetNode(fNSize, ...);
+  }
+  fLambda = n->GetNpar();
+  return kTRUE;
 }
 
 //_________________________________________________________________
@@ -73,23 +100,23 @@ TKDInterpolatorBase::~TKDInterpolatorBase()
     for(int idim=0; idim<fNSize; idim++) delete [] fRefPoints[idim] ;
     delete [] fRefPoints;
   }
-  if(fTNodes){ 
-    fTNodes->Delete();
-    delete fTNodes;
+  if(fNodes){ 
+    fNodes->Delete();
+    delete fNodes;
   }
-  delete [] fTNodesDraw;
+  if(fNodesDraw) delete [] fNodesDraw;
 
-  TH2 *h2=0x0;
-  if((h2 = (TH2S*)gROOT->FindObject("hKDnodes"))) delete h2; 
+  TH2 *h2=NULL;
+  if((h2 = (TH2S*)gROOT->FindObject("hKDnodes"))) delete h2;
 }
 
 
 //__________________________________________________________________
 Bool_t TKDInterpolatorBase::GetCOGPoint(Int_t inode, Float_t *&coord, Float_t &val, Float_t &err) const
 {
-  if(inode < 0 || inode > fNTNodes) return kFALSE;
+  if(inode < 0 || inode > GetNTNodes()) return kFALSE;
 
-  TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode];
+  TKDNodeInfo *node = (TKDNodeInfo*)(*fNodes)[inode];
   coord = &(node->Data()[0]);
   val = node->Val()[0];
   err = node->Val()[1];
@@ -99,23 +126,29 @@ Bool_t     TKDInterpolatorBase::GetCOGPoint(Int_t inode, Float_t *&coord, Float_t &v
 //_________________________________________________________________
 TKDNodeInfo* TKDInterpolatorBase::GetNodeInfo(Int_t inode) const
 {
-  if(!fTNodes || inode >= fNTNodes) return 0x0;
-  return (TKDNodeInfo*)(*fTNodes)[inode];
+  if(!fNodes || inode >= GetNTNodes()) return NULL;
+  return (TKDNodeInfo*)(*fNodes)[inode];
+}
+
+//_________________________________________________________________
+Int_t TKDInterpolatorBase::GetNTNodes() const 
+{
+  return fNodes?fNodes->GetEntriesFast():0;
 }
 
 //_________________________________________________________________
 Bool_t TKDInterpolatorBase::GetRange(Int_t ax, Float_t &min, Float_t &max) const
 {
-  if(!fTNodes) return kFALSE;
-  Int_t ndim = ((TKDNodeInfo*)(*fTNodes)[0])->GetDimension();
+  if(!fNodes) return kFALSE;
+  Int_t ndim = ((TKDNodeInfo*)(*fNodes)[0])->GetDimension();
   if(ax<0 || ax>=ndim){
     min=0.; max=0.;
     return kFALSE;
   }
   min=1.e10; max=-1.e10;
   Float_t axmin, axmax;
-  for(Int_t in=fNTNodes; in--; ){ 
-    TKDNodeInfo *node = (TKDNodeInfo*)((*fTNodes)[in]);
+  for(Int_t in=GetNTNodes(); in--; ){ 
+    TKDNodeInfo *node = (TKDNodeInfo*)((*fNodes)[in]);
     node->GetBoundary(ax, axmin, axmax);
     if(axmin<min) min = axmin;
     if(axmax>max) max = axmax;
@@ -136,9 +169,9 @@ void TKDInterpolatorBase::GetStatus(Option_t *opt)
   printf("  Weights: %s\n", UseWeights() ? "YES" : "NO");
   
   if(strcmp(opt, "all") != 0 ) return;
-  printf("fNTNodes %d\n", fNTNodes);        //Number of evaluation data points
-  for(int i=0; i<fNTNodes; i++){
-    TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[i]; 
+  printf("GetNTNodes() %d\n", GetNTNodes());        //Number of evaluation data points
+  for(int i=0; i<GetNTNodes(); i++){
+    TKDNodeInfo *node = (TKDNodeInfo*)(*fNodes)[i]; 
     printf("%d ", i); node->Print();
   }
 }
@@ -152,37 +185,40 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub
 //
 // 1. The default method used for interpolation is kCOG.
 // 2. The initial number of neighbors used for the estimation is set to Int(alpha*fLambda) (alpha = 1.5)
-      
+                   
   Float_t pointF[50]; // local Float_t conversion for "point"
   for(int idim=0; idim<fNSize; idim++) pointF[idim] = (Float_t)point[idim];
   Int_t nodeIndex = GetNodeIndex(pointF);
-  if(nodeIndex<0){
-    Error("TKDInterpolatorBase::Eval()", "Can not retrive node for data point.");
-    result = 0.;
+  if(nodeIndex<0){ 
+    Error("TKDInterpolatorBase::Eval()", "Can not retrive node for data point.");      
+    result = 0.;   
     error = 1.E10;
     return 0.;
   }
-  TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[nodeIndex];
-  if(node->Cov() && !force) return node->CookPDF(point, result, error);
+  TKDNodeInfo *node = (TKDNodeInfo*)(*fNodes)[nodeIndex];
+  if(node->Par() && !force){ 
+    //printf("Node @ %d\n", nodeIndex); node->Print("a");
+    return node->CookPDF(point, result, error);
+  }
 
   // Allocate memory
   if(!fBuffer) fBuffer = new Double_t[2*fLambda];
   if(!fKDhelper){ 
     fRefPoints = new Float_t*[fNSize];
     for(int id=0; id<fNSize; id++){
-      fRefPoints[id] = new Float_t[fNTNodes];
-      for(int in=0; in<fNTNodes; in++) fRefPoints[id][in] = ((TKDNodeInfo*)(*fTNodes)[in])->Data()[id];
+      fRefPoints[id] = new Float_t[GetNTNodes()];
+      for(int in=0; in<GetNTNodes(); in++) fRefPoints[id][in] = ((TKDNodeInfo*)(*fNodes)[in])->Data()[id];
     }
-    Info("TKDInterpolatorBase::Eval()", Form("Build TKDTree(%d, %d, %d)", fNTNodes, fNSize, kNhelper));
-    fKDhelper = new TKDTreeIF(fNTNodes, fNSize, kNhelper, fRefPoints);
+    Info("TKDInterpolatorBase::Eval()", Form("Build TKDTree(%d, %d, %d)", GetNTNodes(), fNSize, kNhelper));
+    fKDhelper = new TKDTreeIF(GetNTNodes(), fNSize, kNhelper, fRefPoints);
     fKDhelper->Build();
     fKDhelper->MakeBoundariesExact();
   }
   if(!fFitter) fFitter = new TLinearFitter(fLambda, Form("hyp%d", fLambda-1));
   
   // generate parabolic for nD
-  //Float_t alpha = Float_t(2*lambda + 1) / fNTNodes; // the bandwidth or smoothing parameter
-  //Int_t npoints = Int_t(alpha * fNTNodes);
+  //Float_t alpha = Float_t(2*lambda + 1) / GetNTNodes(); // the bandwidth or smoothing parameter
+  //Int_t npoints = Int_t(alpha * GetNTNodes());
   //printf("Params : %d NPoints %d\n", lambda, npoints);
   // prepare workers
 
@@ -206,9 +242,9 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub
       error = 1.E10;
       return 0.;
     } else npoints = npoints_new;
-    if(npoints > fNTNodes){
-      Warning("TKDInterpolatorBase::Eval()", Form("The number of interpolation points requested (%d) exceeds number of PDF values (%d). Downscale.", npoints, fNTNodes));
-      npoints = fNTNodes;
+    if(npoints > GetNTNodes()){
+      Warning("TKDInterpolatorBase::Eval()", Form("The number of interpolation points requested (%d) exceeds number of PDF values (%d). Downscale.", npoints, GetNTNodes()));
+      npoints = GetNTNodes();
       kDOWN = kTRUE;
     }
 
@@ -218,9 +254,9 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub
 
     // add points to fitter
     fFitter->ClearPoints();
-    TKDNodeInfo *tnode = 0x0;
+    TKDNodeInfo *tnode = NULL;
     for(int in=0; in<npoints; in++){
-      tnode = (TKDNodeInfo*)(*fTNodes)[index[in]];
+      tnode = (TKDNodeInfo*)(*fNodes)[index[in]];
       //tnode->Print();
       if(UseCOG()){ // COG
         Float_t *p = &(tnode->Data()[0]);
@@ -264,7 +300,7 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub
   Double_t chi2 = fFitter->GetChisquare()/(npoints - 4 - fLambda);
 
   // store results
-  if(HasStore()) node->Store(par, cov);
+  node->Store(&par, HasStore()?&cov:NULL);
     
   // Build df/dpi|x values
   Double_t *fdfdp = &fBuffer[fLambda];
@@ -282,7 +318,6 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub
     for(int j=0; j<fLambda; j++) error += fdfdp[i]*fdfdp[j]*cov(i,j);
   }    
   error = TMath::Sqrt(error);
-
   return chi2;
 }
 
@@ -297,7 +332,7 @@ void TKDInterpolatorBase::DrawProjection(UInt_t ax1, UInt_t ax2)
   Float_t ax1min, ax1max, ax2min, ax2max;
   GetRange(ax1, ax1min, ax1max);
   GetRange(ax2, ax2min, ax2max);
-  TH2 *h2 = 0x0;
+  TH2 *h2 = NULL;
   if(!(h2 = (TH2S*)gROOT->FindObject("hKDnodes"))){
     h2 = new TH2S("hKDnodes", "", 100, ax1min, ax1max, 100, ax2min, ax2max);
   }
@@ -308,11 +343,11 @@ void TKDInterpolatorBase::DrawProjection(UInt_t ax1, UInt_t ax2)
   h2->Draw();
 
 
-  if(!fTNodesDraw) fTNodesDraw = new TKDNodeInfo::TKDNodeDraw[fNTNodes]; 
-  TKDNodeInfo::TKDNodeDraw *box = 0x0;
-  for(Int_t in=fNTNodes; in--; ){ 
-    box = &(fTNodesDraw[in]);
-    box->SetNode((TKDNodeInfo*)((*fTNodes)[in]), fNSize, ax1, ax2);
+  if(!fNodesDraw) fNodesDraw = new TKDNodeInfo::TKDNodeDraw[GetNTNodes()]; 
+  TKDNodeInfo::TKDNodeDraw *box = NULL;
+  for(Int_t in=GetNTNodes(); in--; ){ 
+    box = &(fNodesDraw[in]);
+    box->SetNode((TKDNodeInfo*)((*fNodes)[in]), fNSize, ax1, ax2);
     box->Draw();
   }
 
@@ -328,10 +363,10 @@ void TKDInterpolatorBase::SetAlpha(Float_t a)
     return;
   }
   // check value
-  if(Int_t((a+1.)*fLambda) > fNTNodes){
-    fAlpha = TMath::Max(0.5, Float_t(fNTNodes)/fLambda-1.);
+  if(Int_t((a+1.)*fLambda) > GetNTNodes()){
+    fAlpha = TMath::Max(0.5, Float_t(GetNTNodes())/fLambda-1.);
     Warning("TKDInterpolatorBase::SetAlpha()", Form("Interpolation neighborhood  exceeds number of evaluation points. Downscale alpha to %f", fAlpha));
-    printf("n[%d] nodes[%d]\n", Int_t((fAlpha+1.)*fLambda), fNTNodes);
+    printf("n[%d] nodes[%d]\n", Int_t((fAlpha+1.)*fLambda), GetNTNodes());
     return;
   }
   fAlpha = a;
index bb8c451..f21e08d 100644 (file)
@@ -41,41 +41,39 @@ public:
   TKDInterpolatorBase(Int_t size = 0);
   virtual ~TKDInterpolatorBase();
 
-  Double_t   Eval(const Double_t *point, Double_t &result, Double_t &error, Bool_t force = kFALSE);
+  Bool_t        Bootstrap();
+  Double_t      Eval(const Double_t *point, Double_t &result, Double_t &error, Bool_t force = kFALSE);
   virtual Int_t GetNodeIndex(const Float_t *p) = 0;
-  Float_t    GetAlpha() const {return fAlpha;}
-  Bool_t     GetCOGPoint(Int_t node, Float_t *&coord, Float_t &val, Float_t &error) const;
-  TKDNodeInfo* GetNodeInfo(Int_t inode) const;
-  Int_t      GetNTNodes() const {return fNTNodes;}
-  Bool_t     GetRange(Int_t ax, Float_t &min, Float_t &max) const;
-  void       GetStatus(Option_t *opt="");
-
-  Bool_t     HasStore() const {return TESTBIT(fStatus, kSTORE);}
-  Bool_t     UseCOG() const {return TESTBIT(fStatus, kCOG);}
-  Bool_t     UseWeights() const {return TESTBIT(fStatus, kWEIGHTS);}
-
-  void       DrawProjection(UInt_t ax1 = 0, UInt_t ax2 = 1);
-  void       SetAlpha(Float_t a);
-  void       SetCOG(Bool_t on = kTRUE) {on ? SETBIT(fStatus, kCOG) : CLRBIT(fStatus, kCOG);}
-  void       SetStore(Bool_t on = kTRUE) {on ? SETBIT(fStatus, kSTORE) : CLRBIT(fStatus, kSTORE);}
-  void       SetWeights(Bool_t on = kTRUE) {on ? SETBIT(fStatus, kWEIGHTS) : CLRBIT(fStatus, kWEIGHTS);}
+  Float_t       GetAlpha() const {return fAlpha;}
+  Int_t         GetLambda() const {return fLambda;}
+  Int_t         GetSize() const {return fNSize;}
+  Bool_t        GetCOGPoint(Int_t node, Float_t *&coord, Float_t &val, Float_t &error) const;
+  TKDNodeInfo*  GetNodeInfo(Int_t inode) const;
+  Int_t         GetNTNodes() const;
+  Bool_t        GetRange(Int_t ax, Float_t &min, Float_t &max) const;
+  void          GetStatus(Option_t *opt="");
 
-protected:
-  virtual void      Build(Int_t nnodes);
+  Bool_t        HasStore() const {return TESTBIT(fStatus, kSTORE);}
+  Bool_t        UseCOG() const {return TESTBIT(fStatus, kCOG);}
+  Bool_t        UseWeights() const {return TESTBIT(fStatus, kWEIGHTS);}
 
-private:
-  TKDInterpolatorBase(const TKDInterpolatorBase &);
-  TKDInterpolatorBase& operator=(const TKDInterpolatorBase &);
+  void          DrawProjection(UInt_t ax1 = 0, UInt_t ax2 = 1);
+  void          SetAlpha(Float_t a);
+  void          SetCOG(Bool_t on = kTRUE) {on ? SETBIT(fStatus, kCOG) : CLRBIT(fStatus, kCOG);}
+  void          SetStore(Bool_t on = kTRUE) {on ? SETBIT(fStatus, kSTORE) : CLRBIT(fStatus, kSTORE);}
+  void          SetWeights(Bool_t on = kTRUE) {on ? SETBIT(fStatus, kWEIGHTS) : CLRBIT(fStatus, kWEIGHTS);}
 
 protected:
-  Int_t         fNSize;       // data dimension
-  Int_t         fNTNodes;     //Number of evaluation data points
-  TClonesArray  *fTNodes;     //interpolation nodes
-  TKDNodeInfo::TKDNodeDraw  *fTNodesDraw; //graphical representation of interpolation nodes
+  virtual Bool_t    Build(Int_t nnodes);
+
+
+  Int_t         fNSize;       //!data dimension
+  TClonesArray  *fNodes;     //interpolation nodes
+  TKDNodeInfo::TKDNodeDraw  *fNodesDraw; //!graphical representation of interpolation nodes
 
 //private:
   UChar_t       fStatus;      // status of the interpolator
-  UChar_t       fLambda;      // number of parameters in polynom
+  UChar_t       fLambda;      //! number of parameters in polynom
   Short_t                  fDepth;       //! depth of the KD Tree structure used
   Float_t       fAlpha;       // parameter controlling the size of the region to interpolate n = (1+alpha)*lambda
   Float_t       **fRefPoints; //! temporary storage of COG data
@@ -83,7 +81,11 @@ protected:
   TKDTree<Int_t, Float_t> *fKDhelper;      //! kNN finder
   TLinearFitter *fFitter;     //! linear fitter        
 
-  ClassDef(TKDInterpolatorBase, 1)   // data interpolator based on KD tree
+private:
+  TKDInterpolatorBase(const TKDInterpolatorBase &);
+  TKDInterpolatorBase& operator=(const TKDInterpolatorBase &);
+
+  ClassDef(TKDInterpolatorBase, 3)   // data interpolator based on KD tree
 };
 
 
index 85b25eb..915a688 100644 (file)
@@ -19,9 +19,9 @@
 
 #include "TKDNodeInfo.h"
 
-#include "TRandom.h"
 #include "TVectorD.h"
 #include "TMatrixD.h"
+#include "TRandom.h"
 #include "TMath.h"
 
 ClassImp(TKDNodeInfo)
@@ -32,9 +32,11 @@ ClassImp(TKDNodeInfo::TKDNodeDraw)
 TKDNodeInfo::TKDNodeInfo(Int_t dim):
   TObject()
   ,fNDim(3*dim)
-  ,fData(0x0)
-  ,fCov(0x0)
-  ,fPar(0x0)
+  ,fData(NULL)
+  ,fNpar(0)
+  ,fNcov(0)
+  ,fPar(NULL)
+  ,fCov(NULL)
 {
   // Default constructor
   fVal[0] = 0.; fVal[1] = 0.;
@@ -45,9 +47,11 @@ TKDNodeInfo::TKDNodeInfo(Int_t dim):
 TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
   TObject((TObject&) ref)
   ,fNDim(fNDim)
-  ,fData(0x0)
-  ,fCov(0x0)
-  ,fPar(0x0)
+  ,fData(NULL)
+  ,fNpar(0)
+  ,fNcov(0)
+  ,fPar(NULL)
+  ,fCov(NULL)
 {
   // Copy constructor
   Build(fNDim/3);
@@ -55,8 +59,16 @@ TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
   memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
   fVal[0] = ref.fVal[0];
   fVal[1] = ref.fVal[1];
-  if(ref.fCov) (*fCov) = (*ref.fCov);
-  if(ref.fPar) (*fPar) = (*ref.fPar);
+  if(ref.fNpar&&ref.fPar){ 
+    fNpar = ref.fNpar;
+    fPar=new Double_t[fNpar];
+    memcpy(fPar, ref.fPar, fNpar*sizeof(Double_t));
+  }
+  if(ref.fNcov && ref.fCov){ 
+    fNcov = ref.fNcov;
+    fCov=new Double_t[fNcov];
+    memcpy(fCov, ref.fCov, fNcov*sizeof(Double_t));
+  }
 }
 
 
@@ -65,10 +77,8 @@ TKDNodeInfo::~TKDNodeInfo()
 {
   // Destructor
   if(fData) delete [] fData;
-  if(fCov){
-    delete fPar;
-    delete fCov;
-  }
+  if(fPar) delete [] fPar;
+  if(fCov) delete [] fCov;
 }
 
 //_________________________________________________________________
@@ -84,9 +94,16 @@ TKDNodeInfo& TKDNodeInfo::operator=(const TKDNodeInfo & ref)
   memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
   fVal[0] = ref.fVal[0];
   fVal[1] = ref.fVal[1];
-  if(ref.fCov) (*fCov) = (*ref.fCov);
-  if(ref.fPar) (*fPar) = (*ref.fPar);
-  
+  if(ref.fNpar&&ref.fPar){ 
+    fNpar = ref.fNpar;
+    fPar=new Double_t[fNpar];
+    memcpy(fPar, ref.fPar, fNpar*sizeof(Double_t));
+  }
+  if(ref.fNcov && ref.fCov){ 
+    fNcov = ref.fNcov;
+    fCov=new Double_t[fNcov];
+    memcpy(fCov, ref.fCov, fNcov*sizeof(Double_t));
+  }
   return *this;
 }
 
@@ -98,18 +115,22 @@ void TKDNodeInfo::Build(Int_t dim)
 //     Info("Build()", "...");
 
   if(!dim) return;
-  fNDim = 3*dim;  
-  Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
+  fNDim = 3*dim;
   if(fData) delete [] fData;
   fData = new Float_t[fNDim];
-  if(fCov){
-    fCov->ResizeTo(lambda, lambda);
-    fPar->ResizeTo(lambda);
-  }
   return;
 }
 
 //_________________________________________________________________
+void TKDNodeInfo::Bootstrap()
+{
+  if(!fNpar || !fPar) return;
+
+  Int_t ndim = Int_t(.5*(TMath::Sqrt(1.+8.*fNpar)-1.))-1;
+  fNDim = 3*ndim;
+}
+
+//_________________________________________________________________
 void TKDNodeInfo::SetNode(Int_t ndim, Float_t *data, Float_t *pdf)
 {
   Build(ndim);
@@ -119,52 +140,58 @@ void TKDNodeInfo::SetNode(Int_t ndim, Float_t *data, Float_t *pdf)
 
 
 //_________________________________________________________________
-void TKDNodeInfo::Print(const Option_t *) const
+void TKDNodeInfo::Print(const Option_t *opt) const
 {
   // Print the content of the node
-  Int_t dim = fNDim/3;
+  Int_t dim = Int_t(fNDim/3.);
   printf("x[");
-  for(int idim=0; idim<dim; idim++) printf("%f ", fData[idim]);
+  for(int idim=0; idim<dim; idim++) printf("%f ", fData?fData[idim]:0.);
   printf("] f = [%f +- %f]\n", fVal[0], fVal[1]);
 
-  Float_t *bounds = &fData[dim];
-  printf("range[");
-  for(int idim=0; idim<dim; idim++) printf("(%f %f) ", bounds[2*idim], bounds[2*idim+1]);
-  printf("]\n");
-  
-  printf("Fit parameters : ");
-  if(!fCov){
-    printf("Not defined.\n");
-    return;
+  if(fData){
+    Float_t *bounds = &fData[dim];
+    printf("range[");
+    for(int idim=0; idim<dim; idim++) printf("(%f %f) ", bounds[2*idim], bounds[2*idim+1]);
+    printf("]\n");
+  }
+  if(strcmp(opt, "a")!=0) return;
+
+  if(fNpar){ 
+    printf("Fit parameters : \n");
+    for(int ip=0; ip<fNpar; ip++) printf("p%d[%f] ", ip, fPar[ip]);
+    printf("\n");
+  }
+  if(!fNcov) return;
+  for(int ip(0), n(0); ip<fNpar; ip++){
+    for(int jp(ip); jp<fNpar; jp++) printf("c(%d %d)[%f] ", ip, jp, fCov[n++]);
+    printf("\n");
   }
-  
-  //   Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
-  for(int ip=0; ip<3; ip++) printf("p%d[%f] ", ip, (*fPar)(ip));
-  printf("\n");
 }
 
 //_________________________________________________________________
-void TKDNodeInfo::Store(const TVectorD &par, const TMatrixD &cov)
+void TKDNodeInfo::Store(TVectorD const *par, TMatrixD const *cov)
 {
-  // Store the parameters and the covariance in the node
-  if(!fCov){
-    fCov = new TMatrixD(cov.GetNrows(), cov.GetNrows());
-    fPar = new TVectorD(par.GetNrows());
-  }
-  (*fPar) = par;
-  (*fCov) = cov;
+// Store the parameters and the covariance in the node
+
+  if(!fPar){SetNpar(); fPar = new Double_t[fNpar];}
+  for(int ip=0; ip<fNpar; ip++) fPar[ip] = (*par)[ip];
+
+  if(!cov) return;
+  if(!fCov){SetNcov(); fCov = new Double_t[fNcov];}
+  for(int ip(0), np(0); ip<fNpar; ip++)
+    for(int jp=ip; jp<fNpar; jp++) fCov[np++] = (*cov)(ip,jp);
 }
 
 //_________________________________________________________________
-Double_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error)
+Bool_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error) const
 {
 // Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
 
-  Int_t ndim = fNDim/3;
-  if(ndim>10) return 0.; // support only up to 10 dimensions
+  Int_t ndim = Int_t(fNDim/3.);
+  if(ndim>10) return kFALSE; // support only up to 10 dimensions
+  //printf("ndim[%d] npar[%d] ncov[%d]\n", ndim, fNpar, fNcov);
 
-  Int_t lambda = 1 + ndim + (ndim*(ndim+1)>>1);
-  Double_t fdfdp[66];
+  Double_t fdfdp[66]; memset(fdfdp, 0, ndim*sizeof(Double_t));
   Int_t ipar = 0;
   fdfdp[ipar++] = 1.;
   for(int idim=0; idim<ndim; idim++){
@@ -174,15 +201,18 @@ Double_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t
 
   // calculate estimation
   result =0.; error = 0.;
-  for(int i=0; i<lambda; i++){
-    result += fdfdp[i]*(*fPar)(i);
-    for(int j=0; j<lambda; j++) error += fdfdp[i]*fdfdp[j]*(*fCov)(i,j);
+  for(int i=0; i<fNpar; i++) result += fdfdp[i]*fPar[i];
+  if(!fNcov) return kTRUE;
+
+  for(int i(0), n(0); i<fNpar; i++){
+    error += fdfdp[i]*fdfdp[i]*fCov[n++];
+    for(int j(i+1); j<fNpar; j++) error += 2.*fdfdp[i]*fdfdp[j]*fCov[n++];
   }    
   error = TMath::Sqrt(error);
   
   //printf("TKDNodeInfo::CookPDF() : %6.3f +- %6.3f\n", result, error);
 
-  return 0.;
+  return kTRUE;
 }
 
 
@@ -191,7 +221,7 @@ Double_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t
 TKDNodeInfo::TKDNodeDraw::TKDNodeDraw() 
   :TBox()
   ,fCOG()
-  ,fNode(0x0)
+  ,fNode(NULL)
 {
   SetFillStyle(3002);
   SetFillColor(50+Int_t(gRandom->Uniform()*50.));
index 198aaba..304dd4d 100644 (file)
@@ -50,33 +50,42 @@ public:
   TKDNodeInfo(const TKDNodeInfo & ref);
   TKDNodeInfo& operator=(const TKDNodeInfo & ref);
   virtual ~TKDNodeInfo();
-  Double_t      CookPDF(const Double_t *point, Double_t &result, Double_t &error);
+  Bool_t        CookPDF(const Double_t *point, Double_t &result, Double_t &error) const;
   inline void   GetBoundary(Int_t ax, Float_t &min, Float_t &max) const;
   inline void   GetCOG(Float_t* &p) const;
   Int_t         GetDimension() const { return fNDim/3; }
   void          GetPDF(Float_t &pdf, Float_t &error) const {pdf=fVal[0]; error=fVal[1];}
   Int_t         GetSize() const { return fNDim; }
+  Int_t         GetNcov() const { return fNcov; }
+  Int_t         GetNpar() const { return fNpar; }
   inline Bool_t Has(const Float_t *p) const;
   void          Print(const Option_t * = "") const;
-  void          Store(const TVectorD &par, const TMatrixD &cov);
+  void          Store(TVectorD const *par, TMatrixD const *cov=NULL);
 
-  TMatrixD*     Cov() const  { return fCov; }
-  TVectorD*     Par() const  { return fPar; }
+  Double_t*      Cov() const  { return fCov; }
+  Double_t*      Par() const  { return fPar; }
 
   void          SetNode(Int_t ndim, Float_t *data, Float_t *pdf);
 protected:
+  void          Bootstrap();
+  void          Build(Int_t ndim);
+  void          SetNcov() { fNcov=Int_t(.5*fNpar*(fNpar+1.));}
+  void          SetNpar() { Int_t dim=Int_t(fNDim/3); fNpar=Int_t(1 + dim + .5*dim*(dim+1));}
   Float_t*      Data() { return fData;}
   Float_t*      Val()  { return &fVal[0]; }
-  void          Build(Int_t ndim);
+
 
 private:
-  Int_t     fNDim;          // 3 times data dimension
-  Float_t   *fData;         //[fNDim] node's data
-  Float_t   fVal[2];        // measured value for node
-  TMatrixD  *fCov;          // interpolator covariance matrix
-  TVectorD  *fPar;          // interpolator parameters
+  Int_t     fNDim;          //! 3 times data dimension
+  Float_t   *fData;         //![fNDim] node's data
+  Float_t   fVal[2];        //!measured value for node
+  Int_t     fNpar;          //number of parameters
+  Int_t     fNcov;          //number of covarince elements
+  Double_t  *fPar;          //[fNpar] interpolator parameters
+  Double_t  *fCov;          //[fNcov] interpolator covariance matrix
+
 
-  ClassDef(TKDNodeInfo, 1)  // node info for interpolator
+  ClassDef(TKDNodeInfo, 2)  // node info for interpolator
 };
 
 //_____________________________________________________________________
index 9967bd7..c4c49c8 100644 (file)
@@ -80,26 +80,26 @@ TKDPDF::~TKDPDF()
 }
 
 //_________________________________________________________________
-void TKDPDF::Build(Int_t)
+Bool_t TKDPDF::Build(Int_t)
 {
 // Fill interpolator's data array i.e.
 //  - estimation points 
 //  - corresponding PDF values
 
   TKDTreeIF::Build();
-  fNTNodes = fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);/*TKDTreeIF::GetNTNodes();*/
   if(!fBoundaries) MakeBoundaries();
   fLambda = 1 + fNDim + (fNDim*(fNDim+1)>>1);
   //printf("after MakeBoundaries() %d\n", memory());
   
   // allocate interpolation nodes
+  Int_t fNTNodes = fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);/*TKDTreeIF::GetNTNodes();*/
   TKDInterpolatorBase::Build(fNTNodes);
 
-  TKDNodeInfo *node = 0x0;
-  Float_t *bounds = 0x0;
+  TKDNodeInfo *node = NULL;
+  Float_t *bounds = NULL;
   Int_t *indexPoints;
   for(int inode=0, tnode = fNNodes; inode<fNTNodes-1; inode++, tnode++){
-    node = (TKDNodeInfo*)(*fTNodes)[inode];
+    node = (TKDNodeInfo*)(*fNodes)[inode];
     node->Val()[0] =  Float_t(fBucketSize)/fNPoints;
     bounds = GetBoundary(tnode);
     for(int idim=0; idim<fNDim; idim++) node->Val()[0] /= (bounds[2*idim+1] - bounds[2*idim]);
@@ -119,7 +119,7 @@ void TKDPDF::Build(Int_t)
   Int_t counts = fNPoints%fBucketSize;
   counts = counts ? counts : fBucketSize;
   Int_t inode = fNTNodes - 1, tnode = inode + fNNodes;
-  node = (TKDNodeInfo*)(*fTNodes)[inode];
+  node = (TKDNodeInfo*)(*fNodes)[inode];
   node->Val()[0] =  Float_t(counts)/fNPoints;
   bounds = GetBoundary(tnode);
   for(int idim=0; idim<fNDim; idim++){ 
@@ -142,7 +142,9 @@ void TKDPDF::Build(Int_t)
   memcpy(&(node->Data()[fNDim]), bounds, fNDimm*sizeof(Float_t));
 
   delete [] fBoundaries;
-  fBoundaries = 0x0;
+  fBoundaries = NULL;
+
+  return kTRUE;
 }
 
 
@@ -155,7 +157,7 @@ void TKDPDF::DrawNode(Int_t tnode, UInt_t ax1, UInt_t ax2)
 // This function creates some graphical objects
 // but don't delete it. Abusing this function may cause memory leaks !
 
-  if(tnode < 0 || tnode >= fNTNodes){
+  if(tnode < 0 || tnode >= GetNTNodes()){
     Warning("DrawNode()", Form("Terminal node %d outside defined range.", tnode));
     return;
   }
@@ -172,7 +174,7 @@ void TKDPDF::DrawNode(Int_t tnode, UInt_t ax1, UInt_t ax2)
   for(int ip = 0; ip<nPoints; ip++) g->SetPoint(ip, fData[ax1][index[ip]], fData[ax2][index[ip]]);
 
   // draw estimation point
-  TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode];
+  TKDNodeInfo *node = (TKDNodeInfo*)(*fNodes)[inode];
   TMarker *m=new TMarker(node->Data()[ax1], node->Data()[ax2], 20);
   m->SetMarkerColor(2);
   m->SetMarkerSize(1.7);
index 257098e..e384475 100644 (file)
@@ -32,7 +32,7 @@ public:
 private:
   TKDPDF(const TKDPDF &);
   TKDPDF& operator=(const TKDPDF &);
-  void     Build(Int_t ndim = 0);
+  Bool_t    Build(Int_t ndim = 0);
 
           
   ClassDef(TKDPDF, 1)   // data interpolator based on KD tree
index 7a355a4..ecaef9d 100644 (file)
@@ -18,8 +18,3 @@ TKDSpline::TKDSpline(Int_t npoints, Int_t ndim) :
 }
 
 
-//_________________________________________________________________
-void TKDSpline::Build(Int_t)
-{
-}
-
index 37fef2b..0ee636b 100644 (file)
@@ -12,7 +12,7 @@ public:
        TKDSpline(Int_t npoints, Int_t ndim);
 
 private:
-       void            Build(Int_t ndim = 0);
+       Bool_t          Build(Int_t ndim=0){return kTRUE;}
        
 protected: