JETAN/AliKMeansClustering.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 // Implemenatation of the K-Means Clustering Algorithm
  17 // See http://en.wikipedia.org/wiki/K-means_clustering and references therein.
  18 //
  19 // This particular implementation is the so called Soft K-means algorithm.
  20 // It has been modified to work on the cylindrical topology in eta-phi space.
  21 //
  22 // Author: Andreas Morsch (CERN)
  23 // andreas.morsch@cern.ch
  24
  25 #include "AliKMeansClustering.h"
  26 #include <TMath.h>
  27 #include <TRandom.h>
  28 #include <TH1F.h>
  29
  30 ClassImp(AliKMeansClustering)
  31
  32 Double_t AliKMeansClustering::fBeta = 10.;
  33
  34
  35 Int_t AliKMeansClustering::SoftKMeans(Int_t k, Int_t n, const Double_t* x, const Double_t* y, Double_t* mx, Double_t* my , Double_t* rk )
  36 {
  37     //
  38     // The soft K-means algorithm
  39     //
  40     Int_t i,j;
  41     //
  42     // (1) Initialisation of the k means
  43
  44     for (i = 0; i < k; i++) {
  45         mx[i] = 2. * TMath::Pi() * gRandom->Rndm();
  46         my[i] = -1. + 2. * gRandom->Rndm();
  47     }
  48
  49     //
  50     // (2a) The responsibilities
  51     Double_t** r = new Double_t*[n]; // responsibilities
  52     for (j = 0; j < n; j++) {r[j] = new Double_t[k];}
  53     //
  54     // (2b) Normalisation
  55     Double_t* nr   = new Double_t[n];
  56     // (3) Iterations
  57     Int_t nit = 0;
  58
  59     while(1) {
  60         nit++;
  61       //
  62       // Assignment step
  63       //
  64       for (j = 0; j < n; j++) {
  65         nr[j] = 0.;
  66         for (i = 0; i < k; i++) {
  67           r[j][i] = TMath::Exp(- fBeta * d(mx[i], my[i], x[j], y[j]));
  68           nr[j] += r[j][i];
  69         } // mean i
  70       } // data point j
  71
  72       for (j = 0; j < n; j++) {
  73         for (i = 0; i < k; i++) {
  74           r[j][i] /=  nr[j];
  75         } // mean i
  76       } // data point j
  77
  78         //
  79         // Update step
  80       Double_t di = 0;
  81
  82       for (i = 0; i < k; i++) {
  83           Double_t oldx = mx[i];
  84           Double_t oldy = my[i];
  85
  86           mx[i] = x[0];
  87           my[i] = y[0];
  88           rk[i] = r[0][i];
  89
  90         for (j = 1; j < n; j++) {
  91             Double_t xx =  x[j];
  92 //
  93 // Here we have to take into acount the cylinder topology where phi is defined mod 2xpi
  94 // If two coordinates are separated by more than pi in phi one has to be shifted by +/- 2 pi
  95
  96             Double_t dx = mx[i] - x[j];
  97             if (dx >  TMath::Pi()) xx += 2. * TMath::Pi();
  98             if (dx < -TMath::Pi()) xx -= 2. * TMath::Pi();
  99             mx[i] = mx[i] * rk[i] + r[j][i] * xx;
 100             my[i] = my[i] * rk[i] + r[j][i] * y[j];
 101             rk[i] += r[j][i];
 102             mx[i] /= rk[i];
 103             my[i] /= rk[i];
 104             if (mx[i] > 2. * TMath::Pi()) mx[i] -= 2. * TMath::Pi();
 105             if (mx[i] < 0.              ) mx[i] += 2. * TMath::Pi();
 106         } // Data
 107         di += d(mx[i], my[i], oldx, oldy);
 108       } // means
 109         //
 110         // ending condition
 111       if (di < 1.e-8 || nit > 1000) break;
 112     } // while
 113
 114 // Clean-up
 115     delete[] nr;
 116     for (j = 0; j < n; j++) delete[] r[j];
 117     delete[] r;
 118 //
 119     return (nit < 1000);
 120
 121 }
 122
 123 Int_t AliKMeansClustering::SoftKMeans2(Int_t k, Int_t n, Double_t* x, Double_t* y, Double_t* mx, Double_t* my , Double_t* sigma2, Double_t* rk )
 124 {
 125     //
 126     // The soft K-means algorithm
 127     //
 128     Int_t i,j;
 129     //
 130     // (1) Initialisation of the k means using k-means++ recipe
 131     //
 132      OptimalInit(k, n, x, y, mx, my);
 133     //
 134     // (2a) The responsibilities
 135     Double_t** r = new Double_t*[n]; // responsibilities
 136     for (j = 0; j < n; j++) {r[j] = new Double_t[k];}
 137     //
 138     // (2b) Normalisation
 139     Double_t* nr = new Double_t[n];
 140     //
 141     // (2c) Weights
 142     Double_t* pi = new Double_t[k];
 143     //
 144     //
 145     // (2d) Initialise the responsibilties and weights
 146     for (j = 0; j < n; j++) {
 147       nr[j] = 0.;
 148       for (i = 0; i < k; i++) {
 149         r[j][i] = TMath::Exp(- fBeta * d(mx[i], my[i], x[j], y[j]));
 150         nr[j] += r[j][i];
 151       } // mean i
 152     } // data point j
 153
 154     for (i = 0; i < k; i++) {
 155       rk[i]    = 0.;
 156       sigma2[i] = 1./fBeta;
 157
 158       for (j = 0; j < n; j++) {
 159         r[j][i] /=  nr[j];
 160         rk[i] += r[j][i];
 161       } // mean i
 162       pi[i] = rk[i] / Double_t(n);
 163     } // data point j
 164     // (3) Iterations
 165     Int_t nit = 0;
 166
 167     while(1) {
 168         nit++;
 169       //
 170       // Assignment step
 171       //
 172       for (j = 0; j < n; j++) {
 173         nr[j] = 0.;
 174         for (i = 0; i < k; i++) {
 175           r[j][i] = pi[i] * TMath::Exp(- d(mx[i], my[i], x[j], y[j]) / sigma2[i] )
 176             / (2. * sigma2[i] * TMath::Pi() * TMath::Pi());
 177           nr[j] += r[j][i];
 178         } // mean i
 179       } // data point j
 180
 181       for (i = 0; i < k; i++) {
 182         for (j = 0; j < n; j++) {
 183           r[j][i] /=  nr[j];
 184         } // mean i
 185       } // data point j
 186
 187         //
 188         // Update step
 189       Double_t di = 0;
 190
 191       for (i = 0; i < k; i++) {
 192           Double_t oldx = mx[i];
 193           Double_t oldy = my[i];
 194
 195           mx[i] = x[0];
 196           my[i] = y[0];
 197           rk[i] = r[0][i];
 198         for (j = 1; j < n; j++) {
 199             Double_t xx =  x[j];
 200 //
 201 // Here we have to take into acount the cylinder topology where phi is defined mod 2xpi
 202 // If two coordinates are separated by more than pi in phi one has to be shifted by +/- 2 pi
 203
 204             Double_t dx = mx[i] - x[j];
 205             if (dx >  TMath::Pi()) xx += 2. * TMath::Pi();
 206             if (dx < -TMath::Pi()) xx -= 2. * TMath::Pi();
 207             if (r[j][i] > 1.e-15) {
 208               mx[i] = mx[i] * rk[i] + r[j][i] * xx;
 209               my[i] = my[i] * rk[i] + r[j][i] * y[j];
 210               rk[i] += r[j][i];
 211               mx[i] /= rk[i];
 212               my[i] /= rk[i];
 213             }
 214             if (mx[i] > 2. * TMath::Pi()) mx[i] -= 2. * TMath::Pi();
 215             if (mx[i] < 0.              ) mx[i] += 2. * TMath::Pi();
 216         } // Data
 217         di += d(mx[i], my[i], oldx, oldy);
 218
 219       } // means
 220       //
 221       // Sigma
 222       for (i = 0; i < k; i++) {
 223         sigma2[i] = 0.;
 224         for (j = 0; j < n; j++) {
 225           sigma2[i] += r[j][i] * d(mx[i], my[i], x[j], y[j]);
 226         } // Data
 227         sigma2[i] /= rk[i];
 228         if (sigma2[i] < 0.0025) sigma2[i] = 0.0025;
 229       } // Clusters
 230       //
 231       // Fractions
 232       for (i = 0; i < k; i++) pi[i] = rk[i] / Double_t(n);
 233       //
 234 // ending condition
 235       if (di < 1.e-8 || nit > 1000) break;
 236     } // while
 237
 238 // Clean-up
 239     delete[] nr;
 240     delete[] pi;
 241     for (j = 0; j < n; j++) delete[] r[j];
 242     delete[] r;
 243 //
 244     return (nit < 1000);
 245 }
 246
 247 Int_t AliKMeansClustering::SoftKMeans3(Int_t k, Int_t n, Double_t* x, Double_t* y, Double_t* mx, Double_t* my ,
 248                                        Double_t* sigmax2, Double_t* sigmay2, Double_t* rk )
 249 {
 250     //
 251     // The soft K-means algorithm
 252     //
 253     Int_t i,j;
 254     //
 255     // (1) Initialisation of the k means using k-means++ recipe
 256     //
 257      OptimalInit(k, n, x, y, mx, my);
 258     //
 259     // (2a) The responsibilities
 260     Double_t** r = new Double_t*[n]; // responsibilities
 261     for (j = 0; j < n; j++) {r[j] = new Double_t[k];}
 262     //
 263     // (2b) Normalisation
 264     Double_t* nr = new Double_t[n];
 265     //
 266     // (2c) Weights
 267     Double_t* pi = new Double_t[k];
 268     //
 269     //
 270     // (2d) Initialise the responsibilties and weights
 271     for (j = 0; j < n; j++) {
 272       nr[j] = 0.;
 273       for (i = 0; i < k; i++) {
 274
 275         r[j][i] = TMath::Exp(- fBeta * d(mx[i], my[i], x[j], y[j]));
 276         nr[j] += r[j][i];
 277       } // mean i
 278     } // data point j
 279
 280     for (i = 0; i < k; i++) {
 281       rk[i]    = 0.;
 282       sigmax2[i] = 1./fBeta;
 283       sigmay2[i] = 1./fBeta;
 284
 285       for (j = 0; j < n; j++) {
 286         r[j][i] /=  nr[j];
 287         rk[i] += r[j][i];
 288       } // mean i
 289       pi[i] = rk[i] / Double_t(n);
 290     } // data point j
 291     // (3) Iterations
 292     Int_t nit = 0;
 293
 294     while(1) {
 295         nit++;
 296       //
 297       // Assignment step
 298       //
 299       for (j = 0; j < n; j++) {
 300         nr[j] = 0.;
 301         for (i = 0; i < k; i++) {
 302
 303           Double_t dx = TMath::Abs(mx[i]-x[j]);
 304           if (dx > TMath::Pi()) dx = 2. * TMath::Pi() - dx;
 305           Double_t dy = TMath::Abs(my[i]-y[j]);
 306           r[j][i] = pi[i] * TMath::Exp(-0.5 *  (dx * dx / sigmax2[i] + dy * dy / sigmay2[i]))
 307             / (2. * TMath::Sqrt(sigmax2[i] * sigmay2[i]) * TMath::Pi() * TMath::Pi());
 308           nr[j] += r[j][i];
 309         } // mean i
 310       } // data point j
 311
 312       for (i = 0; i < k; i++) {
 313         for (j = 0; j < n; j++) {
 314           r[j][i] /=  nr[j];
 315         } // mean i
 316       } // data point j
 317
 318         //
 319         // Update step
 320       Double_t di = 0;
 321
 322       for (i = 0; i < k; i++) {
 323           Double_t oldx = mx[i];
 324           Double_t oldy = my[i];
 325
 326           mx[i] = x[0];
 327           my[i] = y[0];
 328           rk[i] = r[0][i];
 329         for (j = 1; j < n; j++) {
 330             Double_t xx =  x[j];
 331 //
 332 // Here we have to take into acount the cylinder topology where phi is defined mod 2xpi
 333 // If two coordinates are separated by more than pi in phi one has to be shifted by +/- 2 pi
 334
 335             Double_t dx = mx[i] - x[j];
 336             if (dx >  TMath::Pi()) xx += 2. * TMath::Pi();
 337             if (dx < -TMath::Pi()) xx -= 2. * TMath::Pi();
 338             if (r[j][i] > 1.e-15) {
 339               mx[i] = mx[i] * rk[i] + r[j][i] * xx;
 340               my[i] = my[i] * rk[i] + r[j][i] * y[j];
 341               rk[i] += r[j][i];
 342               mx[i] /= rk[i];
 343               my[i] /= rk[i];
 344             }
 345             if (mx[i] > 2. * TMath::Pi()) mx[i] -= 2. * TMath::Pi();
 346             if (mx[i] < 0.              ) mx[i] += 2. * TMath::Pi();
 347         } // Data
 348         di += d(mx[i], my[i], oldx, oldy);
 349
 350       } // means
 351       //
 352       // Sigma
 353       for (i = 0; i < k; i++) {
 354         sigmax2[i] = 0.;
 355         sigmay2[i] = 0.;
 356
 357         for (j = 0; j < n; j++) {
 358           Double_t dx = TMath::Abs(mx[i]-x[j]);
 359           if (dx > TMath::Pi()) dx = 2. * TMath::Pi() - dx;
 360           Double_t dy = TMath::Abs(my[i]-y[j]);
 361           sigmax2[i] += r[j][i] * dx * dx;
 362           sigmay2[i] += r[j][i] * dy * dy;
 363         } // Data
 364         sigmax2[i] /= rk[i];
 365         sigmay2[i] /= rk[i];
 366         if (sigmax2[i] < 0.0025) sigmax2[i] = 0.0025;
 367         if (sigmay2[i] < 0.0025) sigmay2[i] = 0.0025;
 368       } // Clusters
 369       //
 370       // Fractions
 371       for (i = 0; i < k; i++) pi[i] = rk[i] / Double_t(n);
 372       //
 373 // ending condition
 374       if (di < 1.e-8 || nit > 1000) break;
 375     } // while
 376
 377 // Clean-up
 378     delete[] nr;
 379     delete[] pi;
 380     for (j = 0; j < n; j++) delete[] r[j];
 381     delete[] r;
 382 //
 383     return (nit < 1000);
 384 }
 385
 386 Double_t AliKMeansClustering::d(Double_t mx, Double_t my, Double_t x, Double_t y)
 387 {
 388     //
 389     // Distance definition
 390     // Quasi - Euclidian on the eta-phi cylinder
 391
 392     Double_t dx = TMath::Abs(mx-x);
 393     if (dx > TMath::Pi()) dx = 2. * TMath::Pi() - dx;
 394
 395     return (0.5*(dx * dx + (my - y) * (my - y)));
 396 }
 397
 398
 399
 400 void AliKMeansClustering::OptimalInit(Int_t k, Int_t n, const Double_t* x, const Double_t* y, Double_t* mx, Double_t* my)
 401 {
 402   //
 403   // Optimal initialisation using the k-means++ algorithm
 404   // http://en.wikipedia.org/wiki/K-means%2B%2B
 405   //
 406   // k-means++ is an algorithm for choosing the initial values for k-means clustering in statistics and machine learning.
 407   // It was proposed in 2007 by David Arthur and Sergei Vassilvitskii as an approximation algorithm for the NP-hard k-means problem---
 408   // a way of avoiding the sometimes poor clusterings found by the standard k-means algorithm.
 409   //
 410   //
 411   TH1F d2("d2", "", n, -0.5, Float_t(n)-0.5);
 412   d2.Reset();
 413
 414   // (1) Chose first center as a random point among the input data.
 415   Int_t ir = Int_t(Float_t(n) * gRandom->Rndm());
 416   mx[0] = x[ir];
 417   my[0] = y[ir];
 418
 419   // (2) Iterate
 420   Int_t icl = 1;
 421   while(icl < k)
 422     {
 423       // find min distance to existing clusters
 424       for (Int_t j = 0; j < n; j++) {
 425         Double_t dmin = 1.e10;
 426         for (Int_t i = 0; i < icl; i++) {
 427           Double_t dij = d(mx[i], my[i], x[j], y[j]);
 428           if (dij < dmin) dmin = dij;
 429         } // clusters
 430         d2.Fill(Float_t(j), dmin);
 431       } // data points
 432       // select a new cluster from data points with probability ~d2
 433       ir = Int_t(d2.GetRandom() + 0.5);
 434       mx[icl] = x[ir];
 435       my[icl] = y[ir];
 436       icl++;
 437     } // icl
 438 }
 439
 440
 441 ClassImp(AliKMeansResult)
 442
 443
 444
 445 AliKMeansResult::AliKMeansResult(Int_t k):
 446     TObject(),
 447     fK(k),
 448     fMx    (new Double_t[k]),
 449     fMy    (new Double_t[k]),
 450     fSigma2(new Double_t[k]),
 451     fRk    (new Double_t[k]),
 452     fTarget(new Double_t[k]),
 453     fInd   (new Int_t[k])
 454 {
 455 // Constructor
 456 }
 457
 458 AliKMeansResult::AliKMeansResult(const AliKMeansResult &res):
 459   TObject(res),
 460   fK(res.GetK()),
 461   fMx(new Double_t[res.GetK()]),
 462   fMy(new Double_t[res.GetK()]),
 463   fSigma2(new Double_t[res.GetK()]),
 464   fRk(new Double_t[res.GetK()]),
 465   fTarget(new Double_t[res.GetK()]),
 466   fInd(new Int_t[res.GetK()])
 467 {
 468   // Copy constructor
 469   for (Int_t i = 0; i <fK; i++) {
 470     fMx[i]     = (res.GetMx())    [i];
 471     fMy[i]     = (res.GetMy())    [i];
 472     fSigma2[i] = (res.GetSigma2())[i];
 473     fRk[i]     = (res.GetRk())    [i];
 474     fTarget[i] = (res.GetTarget())[i];
 475     fInd[i]    = (res.GetInd())   [i];
 476   }
 477 }
 478
 479 AliKMeansResult& AliKMeansResult::operator=(const AliKMeansResult& res)
 480 {
 481   //
 482   // Assignment operator
 483   if (this != &res) {
 484     TObject::operator=(res);
 485     if (fK != res.fK) {
 486         delete [] fMx;
 487         delete [] fMy;
 488         delete [] fSigma2;
 489         delete [] fRk;
 490         delete [] fTarget;
 491         delete [] fInd;
 492         fK = res.fK;
 493         fMx     = new Double_t[fK];
 494         fMy     = new Double_t[fK];
 495         fSigma2 = new Double_t[fK];
 496         fRk     = new Double_t[fK];
 497         fTarget = new Double_t[fK];
 498         fInd    = new    Int_t[fK];
 499     }
 500
 501     fK = res.fK;
 502     memcpy(fMx,     res.fMx,     fK*sizeof(Double_t));
 503     memcpy(fMy,     res.fMy,     fK*sizeof(Double_t));
 504     memcpy(fSigma2, res.fSigma2, fK*sizeof(Double_t));
 505     memcpy(fRk,     res.fRk,     fK*sizeof(Double_t));
 506     memcpy(fTarget, res.fTarget, fK*sizeof(Double_t));
 507     memcpy(fInd,    res.fInd,    fK*sizeof(Int_t));
 508   }
 509   return *this;
 510 }
 511
 512
 513 AliKMeansResult::~AliKMeansResult()
 514 {
 515 // Destructor
 516     delete[] fMx;
 517     delete[] fMy;
 518     delete[] fSigma2;
 519     delete[] fRk;
 520     delete[] fInd;
 521     delete[] fTarget;
 522 }
 523
 524 void AliKMeansResult::Sort()
 525 {
 526   // Build target array and sort
 527   // Sort clusters
 528   for (Int_t i = 0; i < fK; i++) {
 529     if (fRk[i] > 2.9) {
 530       fTarget[i] = fRk[i] / fSigma2[i];
 531     }
 532     else fTarget[i] = 0.;
 533   }
 534
 535   TMath::Sort(fK, fTarget, fInd);
 536 }
 537
 538 void AliKMeansResult::Sort(Int_t n, const Double_t* x, const Double_t* y)
 539 {
 540   // Build target array and sort
 541   for (Int_t i = 0; i < fK; i++)
 542     {
 543       Int_t nc = 0;
 544       for (Int_t j = 0; j < n; j++)
 545         {
 546           if (2. * AliKMeansClustering::d(fMx[i], fMy[i], x[j], y[j])  <  2.28 * fSigma2[i]) nc++;
 547         }
 548
 549       if (nc > 2) {
 550         fTarget[i] = Double_t(nc) / (2.28 * fSigma2[i]);
 551       } else {
 552         fTarget[i] = 0.;
 553       }
 554     }
 555
 556   TMath::Sort(fK, fTarget, fInd);
 557 }
 558
 559 void AliKMeansResult::CopyResults(const AliKMeansResult* res)
 560 {
 561   fK = res->GetK();
 562   for (Int_t i = 0; i <fK; i++) {
 563     fMx[i]     = (res->GetMx())    [i];
 564     fMy[i]     = (res->GetMy())    [i];
 565     fSigma2[i] = (res->GetSigma2())[i];
 566     fRk[i]     = (res->GetRk())    [i];
 567     fTarget[i] = (res->GetTarget())[i];
 568     fInd[i]    = (res->GetInd())   [i];
 569   }
 570 }