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, Double_t* x, 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     Bool_t rmovalStep = kFALSE;
 167
 168     while(1) {
 169         nit++;
 170       //
 171       // Assignment step
 172       //
 173       for (j = 0; j < n; j++) {
 174         nr[j] = 0.;
 175         for (i = 0; i < k; i++) {
 176           r[j][i] = pi[i] * TMath::Exp(- d(mx[i], my[i], x[j], y[j]) / sigma2[i] )
 177             / (2. * sigma2[i] * TMath::Pi() * TMath::Pi());
 178           nr[j] += r[j][i];
 179         } // mean i
 180       } // data point j
 181
 182       for (i = 0; i < k; i++) {
 183         for (j = 0; j < n; j++) {
 184           r[j][i] /=  nr[j];
 185         } // mean i
 186       } // data point j
 187
 188         //
 189         // Update step
 190       Double_t di = 0;
 191
 192       for (i = 0; i < k; i++) {
 193           Double_t oldx = mx[i];
 194           Double_t oldy = my[i];
 195
 196           mx[i] = x[0];
 197           my[i] = y[0];
 198           rk[i] = r[0][i];
 199         for (j = 1; j < n; j++) {
 200             Double_t xx =  x[j];
 201 //
 202 // Here we have to take into acount the cylinder topology where phi is defined mod 2xpi
 203 // If two coordinates are separated by more than pi in phi one has to be shifted by +/- 2 pi
 204
 205             Double_t dx = mx[i] - x[j];
 206             if (dx >  TMath::Pi()) xx += 2. * TMath::Pi();
 207             if (dx < -TMath::Pi()) xx -= 2. * TMath::Pi();
 208             if (r[j][i] > 1.e-15) {
 209               mx[i] = mx[i] * rk[i] + r[j][i] * xx;
 210               my[i] = my[i] * rk[i] + r[j][i] * y[j];
 211               rk[i] += r[j][i];
 212               mx[i] /= rk[i];
 213               my[i] /= rk[i];
 214             }
 215             if (mx[i] > 2. * TMath::Pi()) mx[i] -= 2. * TMath::Pi();
 216             if (mx[i] < 0.              ) mx[i] += 2. * TMath::Pi();
 217         } // Data
 218         di += d(mx[i], my[i], oldx, oldy);
 219
 220       } // means
 221       //
 222       // Sigma
 223       for (i = 0; i < k; i++) {
 224         sigma2[i] = 0.;
 225         for (j = 1; j < n; j++) {
 226           sigma2[i] += r[j][i] * d(mx[i], my[i], x[j], y[j]);
 227         } // Data
 228         sigma2[i] /= rk[i];
 229         if (sigma2[i] < 0.0025) sigma2[i] = 0.0025;
 230       } // Clusters
 231       //
 232       // Fractions
 233       for (i = 0; i < k; i++) pi[i] = rk[i] / Double_t(n);
 234       //
 235 // ending condition
 236       if (di < 1.e-8 || nit > 1000) break;
 237     } // while
 238
 239 // Clean-up
 240     delete[] nr;
 241     delete[] pi;
 242     for (j = 0; j < n; j++) delete[] r[j];
 243     delete[] r;
 244 //
 245     return (nit < 1000);
 246 }
 247
 248 Int_t AliKMeansClustering::SoftKMeans3(Int_t k, Int_t n, Double_t* x, Double_t* y, Double_t* mx, Double_t* my ,
 249                                        Double_t* sigmax2, Double_t* sigmay2, Double_t* rk )
 250 {
 251     //
 252     // The soft K-means algorithm
 253     //
 254     Int_t i,j;
 255     //
 256     // (1) Initialisation of the k means using k-means++ recipe
 257     //
 258      OptimalInit(k, n, x, y, mx, my);
 259     //
 260     // (2a) The responsibilities
 261     Double_t** r = new Double_t*[n]; // responsibilities
 262     for (j = 0; j < n; j++) {r[j] = new Double_t[k];}
 263     //
 264     // (2b) Normalisation
 265     Double_t* nr = new Double_t[n];
 266     //
 267     // (2c) Weights
 268     Double_t* pi = new Double_t[k];
 269     //
 270     //
 271     // (2d) Initialise the responsibilties and weights
 272     for (j = 0; j < n; j++) {
 273       nr[j] = 0.;
 274       for (i = 0; i < k; i++) {
 275
 276         r[j][i] = TMath::Exp(- fBeta * d(mx[i], my[i], x[j], y[j]));
 277         nr[j] += r[j][i];
 278       } // mean i
 279     } // data point j
 280
 281     for (i = 0; i < k; i++) {
 282       rk[i]    = 0.;
 283       sigmax2[i] = 1./fBeta;
 284       sigmay2[i] = 1./fBeta;
 285
 286       for (j = 0; j < n; j++) {
 287         r[j][i] /=  nr[j];
 288         rk[i] += r[j][i];
 289       } // mean i
 290       pi[i] = rk[i] / Double_t(n);
 291     } // data point j
 292     // (3) Iterations
 293     Int_t nit = 0;
 294     Bool_t rmovalStep = kFALSE;
 295
 296     while(1) {
 297         nit++;
 298       //
 299       // Assignment step
 300       //
 301       for (j = 0; j < n; j++) {
 302         nr[j] = 0.;
 303         for (i = 0; i < k; i++) {
 304
 305           Double_t dx = TMath::Abs(mx[i]-x[j]);
 306           if (dx > TMath::Pi()) dx = 2. * TMath::Pi() - dx;
 307           Double_t dy = TMath::Abs(my[i]-y[j]);
 308           r[j][i] = pi[i] * TMath::Exp(-0.5 *  (dx * dx / sigmax2[i] + dy * dy / sigmay2[i]))
 309             / (2. * TMath::Sqrt(sigmax2[i] * sigmay2[i]) * TMath::Pi() * TMath::Pi());
 310           nr[j] += r[j][i];
 311         } // mean i
 312       } // data point j
 313
 314       for (i = 0; i < k; i++) {
 315         for (j = 0; j < n; j++) {
 316           r[j][i] /=  nr[j];
 317         } // mean i
 318       } // data point j
 319
 320         //
 321         // Update step
 322       Double_t di = 0;
 323
 324       for (i = 0; i < k; i++) {
 325           Double_t oldx = mx[i];
 326           Double_t oldy = my[i];
 327
 328           mx[i] = x[0];
 329           my[i] = y[0];
 330           rk[i] = r[0][i];
 331         for (j = 1; j < n; j++) {
 332             Double_t xx =  x[j];
 333 //
 334 // Here we have to take into acount the cylinder topology where phi is defined mod 2xpi
 335 // If two coordinates are separated by more than pi in phi one has to be shifted by +/- 2 pi
 336
 337             Double_t dx = mx[i] - x[j];
 338             if (dx >  TMath::Pi()) xx += 2. * TMath::Pi();
 339             if (dx < -TMath::Pi()) xx -= 2. * TMath::Pi();
 340             if (r[j][i] > 1.e-15) {
 341               mx[i] = mx[i] * rk[i] + r[j][i] * xx;
 342               my[i] = my[i] * rk[i] + r[j][i] * y[j];
 343               rk[i] += r[j][i];
 344               mx[i] /= rk[i];
 345               my[i] /= rk[i];
 346             }
 347             if (mx[i] > 2. * TMath::Pi()) mx[i] -= 2. * TMath::Pi();
 348             if (mx[i] < 0.              ) mx[i] += 2. * TMath::Pi();
 349         } // Data
 350         di += d(mx[i], my[i], oldx, oldy);
 351
 352       } // means
 353       //
 354       // Sigma
 355       for (i = 0; i < k; i++) {
 356         sigmax2[i] = 0.;
 357         sigmay2[i] = 0.;
 358
 359         for (j = 1; j < n; j++) {
 360           Double_t dx = TMath::Abs(mx[i]-x[j]);
 361           if (dx > TMath::Pi()) dx = 2. * TMath::Pi() - dx;
 362           Double_t dy = TMath::Abs(my[i]-y[j]);
 363           sigmax2[i] += r[j][i] * dx * dx;
 364           sigmay2[i] += r[j][i] * dy * dy;
 365         } // Data
 366         sigmax2[i] /= rk[i];
 367         sigmay2[i] /= rk[i];
 368         if (sigmax2[i] < 0.0025) sigmax2[i] = 0.0025;
 369         if (sigmay2[i] < 0.0025) sigmay2[i] = 0.0025;
 370       } // Clusters
 371       //
 372       // Fractions
 373       for (i = 0; i < k; i++) pi[i] = rk[i] / Double_t(n);
 374       //
 375 // ending condition
 376       if (di < 1.e-8 || nit > 1000) break;
 377     } // while
 378
 379 // Clean-up
 380     delete[] nr;
 381     delete[] pi;
 382     for (j = 0; j < n; j++) delete[] r[j];
 383     delete[] r;
 384 //
 385     return (nit < 1000);
 386 }
 387
 388 Double_t AliKMeansClustering::d(Double_t mx, Double_t my, Double_t x, Double_t y)
 389 {
 390     //
 391     // Distance definition
 392     // Quasi - Euclidian on the eta-phi cylinder
 393
 394     Double_t dx = TMath::Abs(mx-x);
 395     if (dx > TMath::Pi()) dx = 2. * TMath::Pi() - dx;
 396
 397     return (0.5*(dx * dx + (my - y) * (my - y)));
 398 }
 399
 400
 401
 402 void AliKMeansClustering::OptimalInit(Int_t k, Int_t n, Double_t* x, Double_t* y, Double_t* mx, Double_t* my)
 403 {
 404   //
 405   // Optimal initialisation using the k-means++ algorithm
 406   // http://en.wikipedia.org/wiki/K-means%2B%2B
 407   //
 408   // k-means++ is an algorithm for choosing the initial values for k-means clustering in statistics and machine learning.
 409   // It was proposed in 2007 by David Arthur and Sergei Vassilvitskii as an approximation algorithm for the NP-hard k-means problem---
 410   // a way of avoiding the sometimes poor clusterings found by the standard k-means algorithm.
 411   //
 412   //
 413   TH1F d2("d2", "", n, -0.5, Float_t(n)-0.5);
 414   d2.Reset();
 415
 416   // (1) Chose first center as a random point among the input data.
 417   Int_t ir = Int_t(Float_t(n) * gRandom->Rndm());
 418   mx[0] = x[ir];
 419   my[0] = y[ir];
 420
 421   // (2) Iterate
 422   Int_t icl = 1;
 423   while(icl < k)
 424     {
 425       // find min distance to existing clusters
 426       for (Int_t j = 0; j < n; j++) {
 427         Double_t dmin = 1.e10;
 428         for (Int_t i = 0; i < icl; i++) {
 429           Double_t dij = d(mx[i], my[i], x[j], y[j]);
 430           if (dij < dmin) dmin = dij;
 431         } // clusters
 432         d2.Fill(Float_t(j), dmin);
 433       } // data points
 434       // select a new cluster from data points with probability ~d2
 435       ir = Int_t(d2.GetRandom() + 0.5);
 436       mx[icl] = x[ir];
 437       my[icl] = y[ir];
 438       icl++;
 439     } // icl
 440 }
 441
 442
 443 Int_t AliKMeansClustering::NTwoSigma(Int_t k, Int_t n, Double_t* x, Double_t* y, Double_t* mx, Double_t* my,
 444                                            Double_t* sigma2x, Double_t* sigma2y)
 445 {
 446   //
 447   // Counting the number of data points within 2sigma of a cluster
 448   //
 449   Int_t n2sigma = 0;
 450   for (Int_t j = 0; j < n; j++) {
 451     Double_t nassign = 0;
 452     for (Int_t i = 0; i < k; i++) {
 453       Double_t dx = TMath::Abs(mx[i]-x[j]);
 454       if (dx > TMath::Pi()) dx = 2. * TMath::Pi() - dx;
 455       Double_t dy = TMath::Abs(my[i]-y[j]);
 456       if ((dx * dx / sigma2x[i] + dy * dy / sigma2y[i]) < 4.) nassign++;
 457     } // clusters
 458     if (nassign > 0) n2sigma++;
 459   } // data points
 460   return (n2sigma);
 461 }
 462
 463