Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TKDEFGT.cxx
Go to the documentation of this file.
1 // @(#)root/gl:$Id$
2 // Author: Timur Pocheptsov 2009
3 /*************************************************************************
4  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
5  * All rights reserved. *
6  * *
7  * For the licensing terms see $ROOTSYS/LICENSE. *
8  * For the list of contributors see $ROOTSYS/README/CREDITS. *
9  *************************************************************************/
10 #include <climits>
11 
12 #include "TError.h"
13 #include "TMath.h"
14 
15 #include "TKDEFGT.h"
16 #include "TGL5D.h"
17 
18 // Kernel density estimator based on Fast Gauss Transform.
19 //
20 // Nice implementation of FGT by Sebastien Paris
21 // can be found here:
22 // http://www.mathworks.com/matlabcentral/fileexchange/17438
23 //
24 // This version is based on his work.
25 
26 ////////////////////////////////////////////////////////////////////////////////
27 ///Constructor.
28 
29 TKDEFGT::TKDEFGT()
30  : fDim(0),
31  fP(0),
32  fK(0),
33  fSigma(1.),
34  fPD(0),
35  fModelValid(kFALSE)
36 {
37 }
38 
39 namespace {
40 
41 UInt_t NChooseK(UInt_t n, UInt_t k);
42 
43 }
44 
45 ////////////////////////////////////////////////////////////////////////////////
46 ///Destructor.
47 
48 TKDEFGT::~TKDEFGT()
49 {
50 }
51 
52 ////////////////////////////////////////////////////////////////////////////////
53 ///Calculate coefficients for FGT.
54 
55 void TKDEFGT::BuildModel(const std::vector<Double_t> &sources, Double_t sigma,
56  UInt_t dim, UInt_t p, UInt_t k)
57 {
58  if (!sources.size()) {
59  Warning("TKDEFGT::BuildModel", "Bad input - zero size vector");
60  return;
61  }
62 
63  if (!dim) {
64  Warning("TKDEFGT::BuildModel", "Number of dimensions is zero");
65  return;
66  }
67 
68  if (!p) {
69  Warning("TKDEFGT::BuildModel", "Order of truncation is zero, 8 will be used");
70  p = 8;
71  }
72 
73  fDim = dim;
74  fP = p;
75  const UInt_t nP = UInt_t(sources.size()) / fDim;
76  fK = !k ? UInt_t(TMath::Sqrt(Double_t(nP))) : k;
77  fSigma = sigma;
78  fPD = NChooseK(fP + fDim - 1, fDim);
79 
80  fWeights.assign(nP, 1.);
81  fXC.assign(fDim * fK, 0.);
82  fA_K.assign(fPD * fK, 0.);
83  fIndxc.assign(fK, 0);
84  fIndx.assign(nP, 0);
85  fXhead.assign(fK, 0);
86  fXboxsz.assign(fK, 0);
87  fDistC.assign(nP, 0.);
88  fC_K.assign(fPD, 0.);
89  fHeads.assign(fDim + 1, 0);
90  fCinds.assign(fPD, 0);
91  fDx.assign(fDim, 0.);
92  fProds.assign(fPD, 0.);
93 
94  Kcenter(sources);
95  Compute_C_k();
96  Compute_A_k(sources);
97 
98  fModelValid = kTRUE;
99 }
100 
101 ////////////////////////////////////////////////////////////////////////////////
102 ///Calculate coefficients for FGT.
103 ///Alternative specialized version for data from TTree.
104 
105 void TKDEFGT::BuildModel(const TGL5DDataSet *sources, Double_t sigma, UInt_t p, UInt_t k)
106 {
107  if (!sources->SelectedSize()) {
108  Warning("TKDEFGT::BuildModel", "Bad input - zero size vector");
109  return;
110  }
111 
112  fDim = 3;
113 
114  if (!p) {
115  Warning("TKDEFGT::BuildModel", "Order of truncation is zero, 8 will be used");
116  p = 8;
117  }
118 
119  fP = p;
120  const UInt_t nP = sources->SelectedSize();
121  fK = !k ? UInt_t(TMath::Sqrt(Double_t(nP))) : k;
122  fSigma = sigma;
123  fPD = NChooseK(fP + fDim - 1, fDim);
124 
125  fWeights.assign(nP, 1.);
126  fXC.assign(fDim * fK, 0.);
127  fA_K.assign(fPD * fK, 0.);
128  fIndxc.assign(fK, 0);
129  fIndx.assign(nP, 0);
130  fXhead.assign(fK, 0);
131  fXboxsz.assign(fK, 0);
132  fDistC.assign(nP, 0.);
133  fC_K.assign(fPD, 0.);
134  fHeads.assign(fDim + 1, 0);
135  fCinds.assign(fPD, 0);
136  fDx.assign(fDim, 0.);
137  fProds.assign(fPD, 0.);
138 
139  Kcenter(sources);
140  Compute_C_k();
141  Compute_A_k(sources);
142 
143  fModelValid = kTRUE;
144 }
145 
146 void BuildModel();
147 
148 namespace {
149 
150 Double_t DDist(const Double_t *x , const Double_t *y , Int_t d);
151 Double_t DDist(Double_t x1, Double_t y1, Double_t z1, Double_t x2, Double_t y2, Double_t z2);
152 
153 UInt_t Idmax(const std::vector<Double_t> &x , UInt_t n);
154 
155 }
156 
157 ////////////////////////////////////////////////////////////////////////////////
158 ///Solve kcenter task.
159 
160 void TKDEFGT::Kcenter(const std::vector<double> &x)
161 {
162  //Randomly pick one node as the first center.
163  const UInt_t nP = UInt_t(x.size()) / fDim;
164 
165  UInt_t *indxc = &fIndxc[0];
166  UInt_t ind = 1;
167  *indxc++ = ind;
168 
169  const Double_t *x_j = &x[0];
170  const Double_t *x_ind = &x[0] + ind * fDim;
171 
172  for (UInt_t j = 0; j < nP; x_j += fDim, ++j) {
173  fDistC[j] = (j == ind) ? 0. : DDist(x_j , x_ind , fDim);
174  fIndx[j] = 0;
175  }
176 
177  for (UInt_t i = 1 ; i < fK; ++i) {
178  ind = Idmax(fDistC, nP);
179  *indxc++ = ind;
180  x_j = &x[0];
181  x_ind = &x[0] + ind * fDim;
182  for (UInt_t j = 0; j < nP; x_j += fDim, ++j) {
183  const Double_t temp = (j==ind) ? 0.0 : DDist(x_j, x_ind, fDim);
184  if (temp < fDistC[j]) {
185  fDistC[j] = temp;
186  fIndx[j] = i;
187  }
188  }
189  }
190 
191  for (UInt_t i = 0, nd = 0 ; i < nP; i++, nd += fDim) {
192  fXboxsz[fIndx[i]]++;
193  Int_t ibase = fIndx[i] * fDim;
194  for (UInt_t j = 0 ; j < fDim; ++j)
195  fXC[j + ibase] += x[j + nd];
196  }
197 
198  for (UInt_t i = 0, ibase = 0; i < fK ; ++i, ibase += fDim) {
199  const Double_t temp = 1. / fXboxsz[i];
200  for (UInt_t j = 0; j < fDim; ++j)
201  fXC[j + ibase] *= temp;
202  }
203 }
204 
205 ////////////////////////////////////////////////////////////////////////////////
206 ///Solve kcenter task. Version for dim == 3 and data from TTree.
207 ///Randomly pick one node as the first center.
208 
209 void TKDEFGT::Kcenter(const TGL5DDataSet *sources)
210 {
211  const UInt_t nP = sources->SelectedSize();
212 
213  UInt_t *indxc = &fIndxc[0];
214  *indxc++ = 1;
215 
216  {
217  //Block to limit the scope of x_ind etc.
218  const Double_t x_ind = sources->V1(1);
219  const Double_t y_ind = sources->V2(1);
220  const Double_t z_ind = sources->V3(1);
221 
222  for (UInt_t j = 0; j < nP; ++j) {
223  const Double_t x_j = sources->V1(j);
224  const Double_t y_j = sources->V2(j);
225  const Double_t z_j = sources->V3(j);
226  fDistC[j] = (j == 1) ? 0. : DDist(x_j, y_j, z_j, x_ind, y_ind, z_ind);
227  fIndx[j] = 0;
228  }
229  //Block to limit the scope of x_ind etc.
230  }
231 
232  for (UInt_t i = 1 ; i < fK; ++i) {
233  const UInt_t ind = Idmax(fDistC, nP);
234  const Double_t x_ind = sources->V1(ind);
235  const Double_t y_ind = sources->V2(ind);
236  const Double_t z_ind = sources->V3(ind);
237 
238  *indxc++ = ind;
239  for (UInt_t j = 0; j < nP; ++j) {
240  const Double_t x_j = sources->V1(j);
241  const Double_t y_j = sources->V2(j);
242  const Double_t z_j = sources->V3(j);
243 
244  const Double_t temp = (j==ind) ? 0.0 : DDist(x_j, y_j, z_j, x_ind, y_ind, z_ind);
245  if (temp < fDistC[j]) {
246  fDistC[j] = temp;
247  fIndx[j] = i;
248  }
249  }
250  }
251 
252  for (UInt_t i = 0, nd = 0 ; i < nP; i++, nd += fDim) {
253  fXboxsz[fIndx[i]]++;
254  UInt_t ibase = fIndx[i] * fDim;
255  fXC[ibase] += sources->V1(i);
256  fXC[ibase + 1] += sources->V2(i);
257  fXC[ibase + 2] += sources->V3(i);
258  }
259 
260  for (UInt_t i = 0, ibase = 0; i < fK ; ++i, ibase += fDim) {
261  const Double_t temp = 1. / fXboxsz[i];
262  for (UInt_t j = 0; j < fDim; ++j)
263  fXC[j + ibase] *= temp;
264  }
265 }
266 
267 ////////////////////////////////////////////////////////////////////////////////
268 ///Coefficients C_K.
269 
270 void TKDEFGT::Compute_C_k()
271 {
272  fHeads[fDim] = UINT_MAX;
273  fCinds[0] = 0;
274  fC_K[0] = 1.0;
275 
276  for (UInt_t k = 1, t = 1, tail = 1; k < fP; ++k, tail = t) {
277  for (UInt_t i = 0; i < fDim; ++i) {
278  const UInt_t head = fHeads[i];
279  fHeads[i] = t;
280  for (UInt_t j = head; j < tail; ++j, ++t) {
281  fCinds[t] = (j < fHeads[i + 1]) ? fCinds[j] + 1 : 1;
282  fC_K[t] = 2.0 * fC_K[j];
283  fC_K[t] /= fCinds[t];
284  }
285  }
286  }
287 }
288 
289 ////////////////////////////////////////////////////////////////////////////////
290 ///Coefficients A_K.
291 
292 void TKDEFGT::Compute_A_k(const std::vector<Double_t> &x)
293 {
294  const Double_t ctesigma = 1. / fSigma;
295  const UInt_t nP = UInt_t(x.size()) / fDim;
296 
297  for (UInt_t n = 0; n < nP; n++) {
298  UInt_t nbase = n * fDim;
299  UInt_t ix2c = fIndx[n];
300  UInt_t ix2cbase = ix2c * fDim;
301  UInt_t ind = ix2c * fPD;
302  Double_t temp = fWeights[n];
303  Double_t sum = 0.0;
304 
305  for (UInt_t i = 0; i < fDim; ++i) {
306  fDx[i] = (x[i + nbase] - fXC[i + ix2cbase]) * ctesigma;
307  sum += fDx[i] * fDx[i];
308  fHeads[i] = 0;
309  }
310 
311  fProds[0] = TMath::Exp(-sum);
312 
313  for (UInt_t k = 1, t = 1, tail = 1; k < fP; ++k, tail = t) {
314  for (UInt_t i = 0; i < fDim; ++i) {
315  const UInt_t head = fHeads[i];
316  fHeads[i] = t;
317  const Double_t temp1 = fDx[i];
318  for (UInt_t j = head; j < tail; ++j, ++t)
319  fProds[t] = temp1 * fProds[j];
320  }
321  }
322 
323  for (UInt_t i = 0 ; i < fPD ; i++)
324  fA_K[i + ind] += temp * fProds[i];
325  }
326 
327  for (UInt_t k = 0; k < fK; ++k) {
328  const UInt_t ind = k * fPD;
329  for (UInt_t i = 0; i < fPD; ++i)
330  fA_K[i + ind] *= fC_K[i];
331  }
332 }
333 
334 ////////////////////////////////////////////////////////////////////////////////
335 ///Coefficients A_K. Special case for TTree and dim == 3.
336 
337 void TKDEFGT::Compute_A_k(const TGL5DDataSet *sources)
338 {
339  const Double_t ctesigma = 1. / fSigma;
340  const UInt_t nP = sources->SelectedSize();
341 
342  for (UInt_t n = 0; n < nP; n++) {
343  UInt_t ix2c = fIndx[n];
344  UInt_t ix2cbase = ix2c * 3;
345  UInt_t ind = ix2c * fPD;
346  Double_t temp = fWeights[n];
347  Double_t sum = 0.0;
348 
349  fDx[0] = (sources->V1(n) - fXC[ix2cbase]) * ctesigma;
350  fDx[1] = (sources->V2(n) - fXC[ix2cbase + 1]) * ctesigma;
351  fDx[2] = (sources->V3(n) - fXC[ix2cbase + 2]) * ctesigma;
352 
353  sum += (fDx[0] * fDx[0] + fDx[1] * fDx[1] + fDx[2] * fDx[2]);
354  fHeads[0] = fHeads[1] = fHeads[2] = 0;
355 
356  fProds[0] = TMath::Exp(-sum);
357 
358  for (UInt_t k = 1, t = 1, tail = 1; k < fP; ++k, tail = t) {
359  for (UInt_t i = 0; i < 3; ++i) {
360  const UInt_t head = fHeads[i];
361  fHeads[i] = t;
362  const Double_t temp1 = fDx[i];
363  for (UInt_t j = head; j < tail; ++j, ++t)
364  fProds[t] = temp1 * fProds[j];
365  }
366  }
367 
368  for (UInt_t i = 0 ; i < fPD ; i++)
369  fA_K[i + ind] += temp * fProds[i];
370  }
371 
372  for (UInt_t k = 0; k < fK; ++k) {
373  const Int_t ind = k * fPD;
374  for (UInt_t i = 0; i < fPD; ++i)
375  fA_K[i + ind] *= fC_K[i];
376  }
377 }
378 
379 namespace {
380 
381 UInt_t InvNChooseK(UInt_t d, UInt_t cnk);
382 
383 }
384 
385 ////////////////////////////////////////////////////////////////////////////////
386 ///Calculate densities.
387 
388 void TKDEFGT::Predict(const std::vector<Double_t> &ts, std::vector<Double_t> &v, Double_t eval)const
389 {
390  if (!fModelValid) {
391  Error("TKDEFGT::Predict", "Call BuildModel first!");
392  return;
393  }
394 
395  if (!ts.size()) {
396  Warning("TKDEFGT::Predict", "Empty targets vector.");
397  return;
398  }
399 
400  v.assign(ts.size() / fDim, 0.);
401 
402  fHeads.assign(fDim + 1, 0);
403  fDx.assign(fDim, 0.);
404  fProds.assign(fPD, 0.);
405 
406  const Double_t ctesigma = 1. / fSigma;
407  const UInt_t p = InvNChooseK(fDim , fPD);
408  const UInt_t nP = UInt_t(ts.size()) / fDim;
409 
410  for (UInt_t m = 0; m < nP; ++m) {
411  Double_t temp = 0.;
412  const UInt_t mbase = m * fDim;
413 
414  for (UInt_t kn = 0 ; kn < fK ; ++kn) {
415  const UInt_t xbase = kn * fDim;
416  const UInt_t ind = kn * fPD;
417  Double_t sum2 = 0.0;
418  for (UInt_t i = 0; i < fDim ; ++i) {
419  fDx[i] = (ts[i + mbase] - fXC[i + xbase]) * ctesigma;
420  sum2 += fDx[i] * fDx[i];
421  fHeads[i] = 0;
422  }
423 
424  if (sum2 > eval) continue; //skip to next kn
425 
426  fProds[0] = TMath::Exp(-sum2);
427 
428  for (UInt_t k = 1, t = 1, tail = 1; k < p; ++k, tail = t) {
429  for (UInt_t i = 0; i < fDim; ++i) {
430  UInt_t head = fHeads[i];
431  fHeads[i] = t;
432  const Double_t temp1 = fDx[i];
433  for (UInt_t j = head; j < tail; ++j, ++t)
434  fProds[t] = temp1 * fProds[j];
435  }
436  }
437 
438  for (UInt_t i = 0 ; i < fPD ; ++i)
439  temp += fA_K[i + ind] * fProds[i];
440  }
441 
442  v[m] = temp;
443  }
444 
445  Double_t dMin = v[0], dMax = dMin;
446  for (UInt_t i = 1; i < nP; ++i) {
447  dMin = TMath::Min(dMin, v[i]);
448  dMax = TMath::Max(dMax, v[i]);
449  }
450 
451  const Double_t dRange = dMax - dMin;
452  for (UInt_t i = 0; i < nP; ++i)
453  v[i] = (v[i] - dMin) / dRange;
454 
455  dMin = v[0], dMax = dMin;
456  for (UInt_t i = 1; i < nP; ++i) {
457  dMin = TMath::Min(dMin, v[i]);
458  dMax = TMath::Max(dMax, v[i]);
459  }
460 }
461 
462 namespace {
463 
464 ////////////////////////////////////////////////////////////////////////////////
465 ///n is always >= k.
466 
467 UInt_t NChooseK(UInt_t n, UInt_t k)
468 {
469  UInt_t n_k = n - k;
470  if (k < n_k) {
471  k = n_k;
472  n_k = n - k;
473  }
474  UInt_t nchsk = 1;
475  for (UInt_t i = 1; i <= n_k; ++i) {
476  nchsk *= ++k;
477  nchsk /= i;
478  }
479 
480  return nchsk;
481 }
482 
483 ////////////////////////////////////////////////////////////////////////////////
484 
485 Double_t DDist(Double_t x1, Double_t y1, Double_t z1, Double_t x2, Double_t y2, Double_t z2)
486 {
487  return (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2);
488 }
489 
490 ////////////////////////////////////////////////////////////////////////////////
491 
492 Double_t DDist(const Double_t *x , const Double_t *y , Int_t d)
493 {
494  Double_t t = 0., s = 0.;
495 
496  for (Int_t i = 0 ; i < d ; i++) {
497  t = (x[i] - y[i]);
498  s += (t * t);
499  }
500 
501  return s;
502 }
503 
504 ////////////////////////////////////////////////////////////////////////////////
505 
506 UInt_t Idmax(const std::vector<Double_t> &x , UInt_t n)
507 {
508  UInt_t k = 0;
509  Double_t t = -1.0;
510  for (UInt_t i = 0; i < n; ++i) {
511  if (t < x[i]) {
512  t = x[i];
513  k = i;
514  }
515  }
516 
517  return k;
518 }
519 
520 ////////////////////////////////////////////////////////////////////////////////
521 
522 UInt_t InvNChooseK(UInt_t d, UInt_t cnk)
523 {
524  UInt_t cted = 1;
525  for (UInt_t i = 2 ; i <= d ; ++i)
526  cted *= i;
527 
528  const UInt_t cte = cnk * cted;
529  UInt_t p = 2, ctep = 2;
530 
531  for (UInt_t i = p + 1; i < p + d; ++i)
532  ctep *= i;
533 
534  while (ctep != cte) {
535  ctep = ((p + d) * ctep) / p;
536  ++p;
537  }
538 
539  return p;
540 }
541 
542 }//anonymous namespace.