Logo ROOT   6.16/01
Reference Guide
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
29TKDEFGT::TKDEFGT()
30 : fDim(0),
31 fP(0),
32 fK(0),
33 fSigma(1.),
34 fPD(0),
35 fModelValid(kFALSE)
36{
37}
38
39namespace {
40
41UInt_t NChooseK(UInt_t n, UInt_t k);
42
43}
44
45////////////////////////////////////////////////////////////////////////////////
46///Destructor.
47
48TKDEFGT::~TKDEFGT()
49{
50}
51
52////////////////////////////////////////////////////////////////////////////////
53///Calculate coefficients for FGT.
54
55void 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
105void 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
146void BuildModel();
147
148namespace {
149
150Double_t DDist(const Double_t *x , const Double_t *y , Int_t d);
151Double_t DDist(Double_t x1, Double_t y1, Double_t z1, Double_t x2, Double_t y2, Double_t z2);
152
153UInt_t Idmax(const std::vector<Double_t> &x , UInt_t n);
154
155}
156
157////////////////////////////////////////////////////////////////////////////////
158///Solve kcenter task.
159
160void 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
209void 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
270void 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
292void 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
337void 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
379namespace {
380
381UInt_t InvNChooseK(UInt_t d, UInt_t cnk);
382
383}
384
385////////////////////////////////////////////////////////////////////////////////
386///Calculate densities.
387
388void 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
462namespace {
463
464////////////////////////////////////////////////////////////////////////////////
465///n is always >= k.
466
467UInt_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
485Double_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
492Double_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
506UInt_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
522UInt_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.
v
SVector< double, 2 > v
Definition: Dict.h:5
d
#define d(i)
Definition: RSha256.hxx:102
x2
static const double x2[5]
Definition: RooGaussKronrodIntegrator1D.cxx:344
x1
static const double x1[5]
Definition: RooGaussKronrodIntegrator1D.cxx:326
Int_t
int Int_t
Definition: RtypesCore.h:41
UInt_t
unsigned int UInt_t
Definition: RtypesCore.h:42
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:88
Double_t
double Double_t
Definition: RtypesCore.h:55
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:87
Error
void Error(const char *location, const char *msgfmt,...)
Warning
void Warning(const char *location, const char *msgfmt,...)
BuildModel
void BuildModel()
TGL5DDataSet::SelectedSize
UInt_t SelectedSize() const
Size of selected sub-range.
Definition: TGL5D.cxx:156
TGL5DDataSet::V1
Double_t V1(UInt_t ind) const
V1 from sub-range, converted to unit cube.
Definition: TGL5D.cxx:164
TGL5DDataSet::V2
Double_t V2(UInt_t ind) const
V2 from sub-range, converted to unit cube.
Definition: TGL5D.cxx:172
TGL5DDataSet::V3
Double_t V3(UInt_t ind) const
V3 from sub-range, converted to unit cube.
Definition: TGL5D.cxx:180
TKDEFGT::TKDEFGT
TKDEFGT()
Constructor.
Definition: TKDEFGT.cxx:29
TKDEFGT::fHeads
std::vector< UInt_t > fHeads
Definition: TKDEFGT.h:37
TKDEFGT::fPD
UInt_t fPD
Definition: TKDEFGT.h:45
TKDEFGT::Compute_C_k
void Compute_C_k()
Coefficients C_K.
Definition: TKDEFGT.cxx:270
TKDEFGT::fIndxc
std::vector< UInt_t > fIndxc
Definition: TKDEFGT.h:28
TKDEFGT::~TKDEFGT
virtual ~TKDEFGT()
Destructor.
Definition: TKDEFGT.cxx:48
TKDEFGT::Predict
void Predict(const std::vector< Double_t > &targets, std::vector< Double_t > &densities, Double_t e) const
Calculate densities.
Definition: TKDEFGT.cxx:388
TKDEFGT::fCinds
std::vector< UInt_t > fCinds
Definition: TKDEFGT.h:35
TKDEFGT::fXhead
std::vector< UInt_t > fXhead
Definition: TKDEFGT.h:31
TKDEFGT::fC_K
std::vector< Double_t > fC_K
Definition: TKDEFGT.h:34
TKDEFGT::fDim
UInt_t fDim
Definition: TKDEFGT.h:41
TKDEFGT::fDistC
std::vector< Double_t > fDistC
Definition: TKDEFGT.h:33
TKDEFGT::Kcenter
void Kcenter(const std::vector< double > &x)
Solve kcenter task.
Definition: TKDEFGT.cxx:160
TKDEFGT::fK
UInt_t fK
Definition: TKDEFGT.h:43
TKDEFGT::fXC
std::vector< Double_t > fXC
Definition: TKDEFGT.h:26
TKDEFGT::fProds
std::vector< Double_t > fProds
Definition: TKDEFGT.h:39
TKDEFGT::fIndx
std::vector< UInt_t > fIndx
Definition: TKDEFGT.h:30
TKDEFGT::fDx
std::vector< Double_t > fDx
Definition: TKDEFGT.h:38
TKDEFGT::Compute_A_k
void Compute_A_k(const std::vector< Double_t > &x)
Coefficients A_K.
Definition: TKDEFGT.cxx:292
TKDEFGT::fXboxsz
std::vector< UInt_t > fXboxsz
Definition: TKDEFGT.h:32
TKDEFGT::fModelValid
Bool_t fModelValid
Definition: TKDEFGT.h:46
TKDEFGT::fA_K
std::vector< Double_t > fA_K
Definition: TKDEFGT.h:29
TKDEFGT::fSigma
Double_t fSigma
Definition: TKDEFGT.h:44
TKDEFGT::fWeights
std::vector< Double_t > fWeights
Definition: TKDEFGT.h:27
TKDEFGT::BuildModel
void BuildModel(const std::vector< Double_t > &sources, Double_t sigma=1., UInt_t dim=3, UInt_t p=8, UInt_t k=0)
Calculate coefficients for FGT.
Definition: TKDEFGT.cxx:55
TKDEFGT::fP
UInt_t fP
Definition: TKDEFGT.h:42
sigma
const Double_t sigma
Definition: h1analysisProxy.h:11
y
Double_t y[n]
Definition: legend1.C:17
x
Double_t x[n]
Definition: legend1.C:17
n
const Int_t n
Definition: legend1.C:16
TGeoUnit::s
static constexpr double s
Definition: TGeoSystemOfUnits.h:162
TMath::Max
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
TMath::Exp
Double_t Exp(Double_t x)
Definition: TMath.h:715
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:679
TMath::Min
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
m
auto * m
Definition: textangle.C:8
sum
static long int sum(long int i)
Definition: Factory.cxx:2258