[u/mrichter/AliRoot.git] / STAT / TKDInterpolatorBase.cxx

#include "TKDInterpolatorBase.h"
#include "TKDNodeInfo.h"
#include "TKDTree.h"

#include "TClonesArray.h"
#include "TLinearFitter.h"
#include "TTree.h"
#include "TH2.h"
#include "TObjArray.h"
#include "TObjString.h"
#include "TBox.h"
#include "TGraph.h"
#include "TMarker.h"
#include "TMath.h"
#include "TVectorD.h"
#include "TMatrixD.h"

ClassImp(TKDInterpolatorBase)

/////////////////////////////////////////////////////////////////////
// Memory setup of protected data members
// fRefPoints : evaluation point of PDF for each terminal node of underlying KD Tree.
// | 1st terminal node (fNDim point coordinates) | 2nd terminal node (fNDim point coordinates) | ...
//
// fRefValues : evaluation value/error of PDF for each terminal node of underlying KD Tree.
// | 1st terminal node (value) | 2nd terminal node (value) | ... | 1st terminal node (error) | 2nd terminal node (error) | ...
//
// status = |0|0|0|0|0|1(tri-cubic weights)|1(STORE)|1 INT(0 COG )|
/////////////////////////////////////////////////////////////////////


//_________________________________________________________________
TKDInterpolatorBase::TKDInterpolatorBase(Int_t dim) :
	fNSize(dim)
	,fNTNodes(0)
	,fTNodes(0x0)
	,fStatus(4)
	,fLambda(1 + dim + (dim*(dim+1)>>1))
	,fDepth(-1)
	,fAlpha(.5)
	,fRefPoints(0x0)
	,fBuffer(0x0)
	,fKDhelper(0x0)
	,fFitter(0x0)
{
// Default constructor. To be used with care since in this case building
// of data structure is completly left to the user responsability.
}

//_________________________________________________________________
void	TKDInterpolatorBase::Build(Int_t n)
{
	// allocate memory for data

	if(fTNodes) delete fTNodes;
	fNTNodes = n;
	fTNodes = new TClonesArray("TKDNodeInfo", fNTNodes);
	for(int in=0; in<fNTNodes; in++) new ((*fTNodes)[in]) TKDNodeInfo(fNSize);
}

//_________________________________________________________________
TKDInterpolatorBase::~TKDInterpolatorBase()
{
	if(fFitter) delete fFitter;
	if(fKDhelper) delete fKDhelper;
	if(fBuffer) delete [] fBuffer;
	
	if(fRefPoints){
		for(int idim=0; idim<fNSize; idim++) delete [] fRefPoints[idim] ;
		delete [] fRefPoints;
	}
	if(fTNodes) delete fTNodes;
}


//__________________________________________________________________
Bool_t	TKDInterpolatorBase::GetCOGPoint(Int_t inode, Float_t *&coord, Float_t &val, Float_t &err) const
{
	if(inode < 0 || inode > fNTNodes) return kFALSE;

	TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode];
	coord = &(node->Data()[0]);
	val = node->Val()[0];
	err = node->Val()[1];
	return kTRUE;
}

//_________________________________________________________________
TKDNodeInfo* TKDInterpolatorBase::GetNodeInfo(Int_t inode) const
{
	if(!fTNodes || inode >= fNTNodes) return 0x0;
	return (TKDNodeInfo*)(*fTNodes)[inode];
}


//__________________________________________________________________
void TKDInterpolatorBase::GetStatus()
{
// Prints the status of the interpolator

	printf("Interpolator Status :\n");
	printf("  Dim    : %d [%d]\n", fNSize, fLambda);
	printf("  Method : %s\n", fStatus&1 ? "INT" : "COG");
	printf("  Store  : %s\n", fStatus&2 ? "YES" : "NO");
	printf("  Weights: %s\n", fStatus&4 ? "YES" : "NO");
	
	printf("fNTNodes %d\n", fNTNodes);        //Number of evaluation data points
	for(int i=0; i<fNTNodes; i++){
		TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[i]; 
		printf("%d ", i); node->Print();
	}
}

//_________________________________________________________________
Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Double_t &error, Bool_t force)
{
// Evaluate PDF for "point". The result is returned in "result" and error in "error". The function returns the chi2 of the fit.
//
// Observations:
//
// 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){
		result = 0.;
		error = 1.E10;
		return 0.;
	}
	TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[nodeIndex];
	if((fStatus&1) && node->Cov() && !force) 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];
		}
		fKDhelper = new TKDTreeIF(fNTNodes, fNSize, 30, fRefPoints);
		fKDhelper->MakeBoundaries();
	}
	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);
	//printf("Params : %d NPoints %d\n", lambda, npoints);
	// prepare workers

	Int_t *index,  // indexes of NN 
	      ipar,    // local looping variable
				npoints = Int_t((1.+fAlpha)*fLambda); // number of data points used for interpolation
	Float_t *dist, // distances of NN
					d,     // NN normalized distance
					w0,    // work
					w;     // tri-cubic weight function

	do{
		// find nearest neighbors
		for(int idim=0; idim<fNSize; idim++) pointF[idim] = (Float_t)point[idim];
		if(!fKDhelper->FindNearestNeighbors(pointF, npoints+1, index, dist)){
			Error("Eval()", Form("Failed retriving %d neighbours for point:", npoints));
			for(int idim=0; idim<fNSize; idim++) printf("%f ", point[idim]);
			printf("\n");
			return -1;
		}
		// add points to fitter
		fFitter->ClearPoints();
		TKDNodeInfo *tnode = 0x0;
		for(int in=0; in<npoints; in++){
			tnode = (TKDNodeInfo*)(*fTNodes)[index[in]];
			//tnode->Print();
			if(fStatus&1){ // INT
				Float_t *bounds = &(tnode->Data()[fNSize]);
				ipar = 0;
				for(int idim=0; idim<fNSize; idim++){
					fBuffer[ipar++] = .5*(bounds[2*idim] + bounds[2*idim+1]);
					fBuffer[ipar++] = (bounds[2*idim]*bounds[2*idim] + bounds[2*idim] * bounds[2*idim+1] + bounds[2*idim+1] * bounds[2*idim+1])/3.;
					for(int jdim=idim+1; jdim<fNSize; jdim++) fBuffer[ipar++] = (bounds[2*idim] + bounds[2*idim+1]) * (bounds[2*jdim] + bounds[2*jdim+1]) * .25;
				}
			} else { // COG
				Float_t *p = &(tnode->Data()[0]);
				ipar = 0;
				for(int idim=0; idim<fNSize; idim++){
					fBuffer[ipar++] = p[idim];
					for(int jdim=idim; jdim<fNSize; jdim++) fBuffer[ipar++] = p[idim]*p[jdim];
				}
			}

			// calculate tri-cubic weighting function
			if(fStatus&4){
				d = dist[in]/ dist[npoints];
				w0 = (1. - d*d*d); w = w0*w0*w0;
			} else w = 1.;
			 
// 			printf("x[");
// 			for(int idim=0; idim<fLambda-1; idim++) printf("%f ", fBuffer[idim]);
// 			printf("]  v[%f +- %f] (%f, %f)\n", tnode->Val()[0], tnode->Val()[1]/w, tnode->Val()[1], w);
			fFitter->AddPoint(fBuffer, tnode->Val()[0], tnode->Val()[1]/w);
		}
		npoints += 4;
	} while(fFitter->Eval());

	// retrive fitter results
	TMatrixD cov(fLambda, fLambda);
	TVectorD par(fLambda);
	fFitter->GetCovarianceMatrix(cov);
	fFitter->GetParameters(par);
	Double_t chi2 = fFitter->GetChisquare()/(npoints - 4 - fLambda);

	// store results
	if(fStatus&2 && fStatus&1) node->Store(par, cov);
		
	// Build df/dpi|x values
	Double_t *fdfdp = &fBuffer[fLambda];
	ipar = 0;
	fdfdp[ipar++] = 1.;
	for(int idim=0; idim<fNSize; idim++){
		fdfdp[ipar++] = point[idim];
		for(int jdim=idim; jdim<fNSize; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
	}

	// calculate estimation
	result =0.; error = 0.;
	for(int i=0; i<fLambda; i++){
		result += fdfdp[i]*par(i);
		for(int j=0; j<fLambda; j++) error += fdfdp[i]*fdfdp[j]*cov(i,j);
	}	
	error = TMath::Sqrt(error);

	return chi2;
}

//_________________________________________________________________
void TKDInterpolatorBase::DrawBins(UInt_t ax1, UInt_t ax2, Float_t ax1min, Float_t ax1max, Float_t ax2min, Float_t ax2max)
{
// Draw nodes structure projected on plane "ax1:ax2". The parameter
// "depth" specifies the bucket size per node. If depth == -1 draw only
// terminal nodes and evaluation points (default -1 i.e. bucket size per node equal bucket size specified by the user)
//
// Observation:
// This function creates the nodes (TBox) array for the specified depth
// but don't delete it. Abusing this function may cause memory leaks !


	TH2 *h2 = new TH2S("hNodes", "", 100, ax1min, ax1max, 100, ax2min, ax2max);
	h2->GetXaxis()->SetTitle(Form("x_{%d}", ax1));
	h2->GetYaxis()->SetTitle(Form("x_{%d}", ax2));
	h2->Draw();
	
	const Float_t kBorder = 0.;//1.E-4;
	TBox *boxArray = new TBox[fNTNodes], *box;
	Float_t *bounds = 0x0;
	for(int inode = 0; inode < fNTNodes; inode++){
		box = &boxArray[inode];
		box->SetFillStyle(3002);
		box->SetFillColor(50+inode/*Int_t(gRandom->Uniform()*50.)*/);
		
		bounds = &(((TKDNodeInfo*)(*fTNodes)[inode])->Data()[fNSize]);
		box->DrawBox(bounds[2*ax1]+kBorder, bounds[2*ax2]+kBorder, bounds[2*ax1+1]-kBorder, bounds[2*ax2+1]-kBorder);
	}

	// Draw reference points
	TGraph *ref = new TGraph(fNTNodes);
	ref->SetMarkerStyle(3);
	ref->SetMarkerSize(.7);
	ref->SetMarkerColor(2);
	for(int inode = 0; inode < fNTNodes; inode++){
		TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode];
		ref->SetPoint(inode, node->Data()[ax1], node->Data()[ax2]);
	}
	ref->Draw("p");
	return;
}

//__________________________________________________________________
void TKDInterpolatorBase::SetInterpolationMethod(Bool_t on)
{
// Set interpolation bit to "on".
	
	if(on) fStatus += fStatus&1 ? 0 : 1;
	else fStatus += fStatus&1 ? -1 : 0;
}


//_________________________________________________________________
void TKDInterpolatorBase::SetStore(Bool_t on)
{
// Set store bit to "on"
	
	if(on) fStatus += fStatus&2 ? 0 : 2;
	else fStatus += fStatus&2 ? -2 : 0;
}

//_________________________________________________________________
void TKDInterpolatorBase::SetWeights(Bool_t on)
{
// Set weights bit to "on"
	
	if(on) fStatus += fStatus&4 ? 0 : 4;
	else fStatus += fStatus&4 ? -4 : 0;
}
Commit	Line	Data
a3408ed3	1	#include "TKDInterpolatorBase.h"
	2	#include "TKDNodeInfo.h"
	3	#include "TKDTree.h"
	4
	5	#include "TClonesArray.h"
	6	#include "TLinearFitter.h"
	7	#include "TTree.h"
	8	#include "TH2.h"
	9	#include "TObjArray.h"
	10	#include "TObjString.h"
	11	#include "TBox.h"
	12	#include "TGraph.h"
	13	#include "TMarker.h"
	14	#include "TMath.h"
	15	#include "TVectorD.h"
	16	#include "TMatrixD.h"
	17
	18	ClassImp(TKDInterpolatorBase)
	19
	20	/////////////////////////////////////////////////////////////////////
	21	// Memory setup of protected data members
	22	// fRefPoints : evaluation point of PDF for each terminal node of underlying KD Tree.
	23	// \| 1st terminal node (fNDim point coordinates) \| 2nd terminal node (fNDim point coordinates) \| ...
	24	//
	25	// fRefValues : evaluation value/error of PDF for each terminal node of underlying KD Tree.
	26	// \| 1st terminal node (value) \| 2nd terminal node (value) \| ... \| 1st terminal node (error) \| 2nd terminal node (error) \| ...
	27	//
	28	// status = \|0\|0\|0\|0\|0\|1(tri-cubic weights)\|1(STORE)\|1 INT(0 COG )\|
	29	/////////////////////////////////////////////////////////////////////
	30
	31
	32	//_________________________________________________________________
1a439134	33	TKDInterpolatorBase::TKDInterpolatorBase(Int_t dim) :
a3408ed3	34	fNSize(dim)
	35	,fNTNodes(0)
	36	,fTNodes(0x0)
	37	,fStatus(4)
	38	,fLambda(1 + dim + (dim*(dim+1)>>1))
	39	,fDepth(-1)
	40	,fAlpha(.5)
	41	,fRefPoints(0x0)
	42	,fBuffer(0x0)
	43	,fKDhelper(0x0)
	44	,fFitter(0x0)
	45	{
	46	// Default constructor. To be used with care since in this case building
	47	// of data structure is completly left to the user responsability.
	48	}
	49
	50	//_________________________________________________________________
1a439134	51	void TKDInterpolatorBase::Build(Int_t n)
a3408ed3	52	{
	53	// allocate memory for data
	54
	55	if(fTNodes) delete fTNodes;
	56	fNTNodes = n;
	57	fTNodes = new TClonesArray("TKDNodeInfo", fNTNodes);
	58	for(int in=0; in<fNTNodes; in++) new ((*fTNodes)[in]) TKDNodeInfo(fNSize);
	59	}
	60
	61	//_________________________________________________________________
	62	TKDInterpolatorBase::~TKDInterpolatorBase()
	63	{
	64	if(fFitter) delete fFitter;
	65	if(fKDhelper) delete fKDhelper;
	66	if(fBuffer) delete [] fBuffer;
	67
	68	if(fRefPoints){
	69	for(int idim=0; idim<fNSize; idim++) delete [] fRefPoints[idim] ;
	70	delete [] fRefPoints;
	71	}
	72	if(fTNodes) delete fTNodes;
	73	}
	74
	75
	76	//__________________________________________________________________
	77	Bool_t TKDInterpolatorBase::GetCOGPoint(Int_t inode, Float_t *&coord, Float_t &val, Float_t &err) const
	78	{
	79	if(inode < 0 \|\| inode > fNTNodes) return kFALSE;
	80
	81	TKDNodeInfo node = (TKDNodeInfo)(*fTNodes)[inode];
1a439134	82	coord = &(node->Data()[0]);
	83	val = node->Val()[0];
	84	err = node->Val()[1];
a3408ed3	85	return kTRUE;
	86	}
	87
	88	//_________________________________________________________________
1a439134	89	TKDNodeInfo* TKDInterpolatorBase::GetNodeInfo(Int_t inode) const
a3408ed3	90	{
	91	if(!fTNodes \|\| inode >= fNTNodes) return 0x0;
	92	return (TKDNodeInfo)(fTNodes)[inode];
	93	}
	94
	95
	96	//__________________________________________________________________
	97	void TKDInterpolatorBase::GetStatus()
	98	{
	99	// Prints the status of the interpolator
	100
	101	printf("Interpolator Status :\n");
	102	printf(" Dim : %d [%d]\n", fNSize, fLambda);
	103	printf(" Method : %s\n", fStatus&1 ? "INT" : "COG");
	104	printf(" Store : %s\n", fStatus&2 ? "YES" : "NO");
	105	printf(" Weights: %s\n", fStatus&4 ? "YES" : "NO");
	106
	107	printf("fNTNodes %d\n", fNTNodes); //Number of evaluation data points
	108	for(int i=0; i<fNTNodes; i++){
	109	TKDNodeInfo node = (TKDNodeInfo)(*fTNodes)[i];
	110	printf("%d ", i); node->Print();
	111	}
	112	}
	113
	114	//_________________________________________________________________
	115	Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Double_t &error, Bool_t force)
	116	{
	117	// Evaluate PDF for "point". The result is returned in "result" and error in "error". The function returns the chi2 of the fit.
	118	//
	119	// Observations:
	120	//
	121	// 1. The default method used for interpolation is kCOG.
	122	// 2. The initial number of neighbors used for the estimation is set to Int(alpha*fLambda) (alpha = 1.5)
	123
	124	Float_t pointF[50]; // local Float_t conversion for "point"
	125	for(int idim=0; idim<fNSize; idim++) pointF[idim] = (Float_t)point[idim];
	126	Int_t nodeIndex = GetNodeIndex(pointF);
	127	if(nodeIndex<0){
	128	result = 0.;
	129	error = 1.E10;
	130	return 0.;
	131	}
	132	TKDNodeInfo node = (TKDNodeInfo)(*fTNodes)[nodeIndex];
1a439134	133	if((fStatus&1) && node->Cov() && !force) return node->CookPDF(point, result, error);
a3408ed3	134
	135	// Allocate memory
	136	if(!fBuffer) fBuffer = new Double_t[2*fLambda];
	137	if(!fKDhelper){
	138	fRefPoints = new Float_t*[fNSize];
	139	for(int id=0; id<fNSize; id++){
	140	fRefPoints[id] = new Float_t[fNTNodes];
1a439134	141	for(int in=0; in<fNTNodes; in++) fRefPoints[id][in] = ((TKDNodeInfo)(fTNodes)[in])->Data()[id];
a3408ed3	142	}
	143	fKDhelper = new TKDTreeIF(fNTNodes, fNSize, 30, fRefPoints);
	144	fKDhelper->MakeBoundaries();
	145	}
	146	if(!fFitter) fFitter = new TLinearFitter(fLambda, Form("hyp%d", fLambda-1));
	147
	148	// generate parabolic for nD
	149	//Float_t alpha = Float_t(2*lambda + 1) / fNTNodes; // the bandwidth or smoothing parameter
	150	//Int_t npoints = Int_t(alpha * fNTNodes);
	151	//printf("Params : %d NPoints %d\n", lambda, npoints);
	152	// prepare workers
	153
	154	Int_t *index, // indexes of NN
	155	ipar, // local looping variable
	156	npoints = Int_t((1.+fAlpha)*fLambda); // number of data points used for interpolation
	157	Float_t *dist, // distances of NN
	158	d, // NN normalized distance
	159	w0, // work
	160	w; // tri-cubic weight function
	161
	162	do{
	163	// find nearest neighbors
	164	for(int idim=0; idim<fNSize; idim++) pointF[idim] = (Float_t)point[idim];
	165	if(!fKDhelper->FindNearestNeighbors(pointF, npoints+1, index, dist)){
	166	Error("Eval()", Form("Failed retriving %d neighbours for point:", npoints));
	167	for(int idim=0; idim<fNSize; idim++) printf("%f ", point[idim]);
	168	printf("\n");
	169	return -1;
	170	}
	171	// add points to fitter
	172	fFitter->ClearPoints();
	173	TKDNodeInfo *tnode = 0x0;
	174	for(int in=0; in<npoints; in++){
	175	tnode = (TKDNodeInfo)(fTNodes)[index[in]];
	176	//tnode->Print();
	177	if(fStatus&1){ // INT
1a439134	178	Float_t *bounds = &(tnode->Data()[fNSize]);
a3408ed3	179	ipar = 0;
	180	for(int idim=0; idim<fNSize; idim++){
	181	fBuffer[ipar++] = .5(bounds[2idim] + bounds[2*idim+1]);
	182	fBuffer[ipar++] = (bounds[2idim]bounds[2idim] + bounds[2idim] * bounds[2idim+1] + bounds[2idim+1] * bounds[2*idim+1])/3.;
	183	for(int jdim=idim+1; jdim<fNSize; jdim++) fBuffer[ipar++] = (bounds[2idim] + bounds[2idim+1]) * (bounds[2jdim] + bounds[2jdim+1]) * .25;
	184	}
	185	} else { // COG
1a439134	186	Float_t *p = &(tnode->Data()[0]);
a3408ed3	187	ipar = 0;
	188	for(int idim=0; idim<fNSize; idim++){
	189	fBuffer[ipar++] = p[idim];
	190	for(int jdim=idim; jdim<fNSize; jdim++) fBuffer[ipar++] = p[idim]*p[jdim];
	191	}
	192	}
	193
	194	// calculate tri-cubic weighting function
	195	if(fStatus&4){
	196	d = dist[in]/ dist[npoints];
	197	w0 = (1. - ddd); w = w0w0w0;
	198	} else w = 1.;
	199
	200	// printf("x[");
	201	// for(int idim=0; idim<fLambda-1; idim++) printf("%f ", fBuffer[idim]);
1a439134	202	// printf("] v[%f +- %f] (%f, %f)\n", tnode->Val()[0], tnode->Val()[1]/w, tnode->Val()[1], w);
1a439134	203	fFitter->AddPoint(fBuffer, tnode->Val()[0], tnode->Val()[1]/w);
a3408ed3	204	}
	205	npoints += 4;
	206	} while(fFitter->Eval());
	207
	208	// retrive fitter results
	209	TMatrixD cov(fLambda, fLambda);
	210	TVectorD par(fLambda);
	211	fFitter->GetCovarianceMatrix(cov);
	212	fFitter->GetParameters(par);
	213	Double_t chi2 = fFitter->GetChisquare()/(npoints - 4 - fLambda);
	214
	215	// store results
	216	if(fStatus&2 && fStatus&1) node->Store(par, cov);
	217
	218	// Build df/dpi\|x values
	219	Double_t *fdfdp = &fBuffer[fLambda];
	220	ipar = 0;
	221	fdfdp[ipar++] = 1.;
	222	for(int idim=0; idim<fNSize; idim++){
	223	fdfdp[ipar++] = point[idim];
	224	for(int jdim=idim; jdim<fNSize; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
	225	}
	226
	227	// calculate estimation
	228	result =0.; error = 0.;
	229	for(int i=0; i<fLambda; i++){
	230	result += fdfdp[i]*par(i);
	231	for(int j=0; j<fLambda; j++) error += fdfdp[i]fdfdp[j]cov(i,j);
	232	}
	233	error = TMath::Sqrt(error);
	234
	235	return chi2;
	236	}
	237
	238	//_________________________________________________________________
	239	void TKDInterpolatorBase::DrawBins(UInt_t ax1, UInt_t ax2, Float_t ax1min, Float_t ax1max, Float_t ax2min, Float_t ax2max)
	240	{
	241	// Draw nodes structure projected on plane "ax1:ax2". The parameter
	242	// "depth" specifies the bucket size per node. If depth == -1 draw only
	243	// terminal nodes and evaluation points (default -1 i.e. bucket size per node equal bucket size specified by the user)
	244	//
	245	// Observation:
	246	// This function creates the nodes (TBox) array for the specified depth
	247	// but don't delete it. Abusing this function may cause memory leaks !
	248
	249
	250
	251	TH2 *h2 = new TH2S("hNodes", "", 100, ax1min, ax1max, 100, ax2min, ax2max);
	252	h2->GetXaxis()->SetTitle(Form("x_{%d}", ax1));
	253	h2->GetYaxis()->SetTitle(Form("x_{%d}", ax2));
	254	h2->Draw();
	255
	256	const Float_t kBorder = 0.;//1.E-4;
	257	TBox boxArray = new TBox[fNTNodes], box;
	258	Float_t *bounds = 0x0;
	259	for(int inode = 0; inode < fNTNodes; inode++){
	260	box = &boxArray[inode];
	261	box->SetFillStyle(3002);
	262	box->SetFillColor(50+inode/Int_t(gRandom->Uniform()50.)*/);
	263
1a439134	264	bounds = &(((TKDNodeInfo)(fTNodes)[inode])->Data()[fNSize]);
a3408ed3	265	box->DrawBox(bounds[2ax1]+kBorder, bounds[2ax2]+kBorder, bounds[2ax1+1]-kBorder, bounds[2ax2+1]-kBorder);
	266	}
	267
	268	// Draw reference points
	269	TGraph *ref = new TGraph(fNTNodes);
	270	ref->SetMarkerStyle(3);
	271	ref->SetMarkerSize(.7);
	272	ref->SetMarkerColor(2);
	273	for(int inode = 0; inode < fNTNodes; inode++){
	274	TKDNodeInfo node = (TKDNodeInfo)(*fTNodes)[inode];
1a439134	275	ref->SetPoint(inode, node->Data()[ax1], node->Data()[ax2]);
a3408ed3	276	}
	277	ref->Draw("p");
	278	return;
	279	}
	280
	281	//__________________________________________________________________
1a439134	282	void TKDInterpolatorBase::SetInterpolationMethod(Bool_t on)
a3408ed3	283	{
	284	// Set interpolation bit to "on".
	285
	286	if(on) fStatus += fStatus&1 ? 0 : 1;
	287	else fStatus += fStatus&1 ? -1 : 0;
	288	}
	289
	290
	291	//_________________________________________________________________
1a439134	292	void TKDInterpolatorBase::SetStore(Bool_t on)
a3408ed3	293	{
	294	// Set store bit to "on"
	295
	296	if(on) fStatus += fStatus&2 ? 0 : 2;
	297	else fStatus += fStatus&2 ? -2 : 0;
	298	}
	299
	300	//_________________________________________________________________
1a439134	301	void TKDInterpolatorBase::SetWeights(Bool_t on)
a3408ed3	302	{
	303	// Set weights bit to "on"
	304
	305	if(on) fStatus += fStatus&4 ? 0 : 4;
	306	else fStatus += fStatus&4 ? -4 : 0;
	307	}