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
  29 ClassImp(AliKMeansClustering)
  30
  31 Double_t AliKMeansClustering::fBeta = 10.;
  32
  33 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 )
  34 {
  35     //
  36     // The soft K-means algorithm
  37     //
  38     Int_t i,j;
  39     //
  40     // (1) Initialisation of the k means
  41
  42     for (i = 0; i < k; i++) {
  43         mx[i] = 2. * TMath::Pi() * gRandom->Rndm();
  44         my[i] = -1. + 2. * gRandom->Rndm();
  45     }
  46
  47     //
  48     // (2a) The responsibilities
  49     Double_t** r = new Double_t*[n]; // responsibilities
  50     for (j = 0; j < n; j++) {r[j] = new Double_t[k];}
  51     //
  52     // (2b) Normalisation
  53     Double_t* nr = new Double_t[n];
  54
  55     // (3) Iterations
  56     Int_t nit = 0;
  57
  58     while(1) {
  59         nit++;
  60       //
  61       // Assignment step
  62       //
  63       for (j = 0; j < n; j++) {
  64         nr[j] = 0.;
  65         for (i = 0; i < k; i++) {
  66           r[j][i] = TMath::Exp(- fBeta * d(mx[i], my[i], x[j], y[j]));
  67           nr[j] += r[j][i];
  68         } // mean i
  69       } // data point j
  70
  71       for (j = 0; j < n; j++) {
  72         for (i = 0; i < k; i++) {
  73           r[j][i] /=  nr[j];
  74         } // mean i
  75       } // data point j
  76
  77         //
  78         // Update step
  79       Double_t di = 0;
  80
  81       for (i = 0; i < k; i++) {
  82           Double_t oldx = mx[i];
  83           Double_t oldy = my[i];
  84
  85           mx[i] = x[0];
  86           my[i] = y[0];
  87           rk[i] = r[0][i];
  88
  89         for (j = 1; j < n; j++) {
  90             Double_t xx =  x[j];
  91 //
  92 // Here we have to take into acount the cylinder topology where phi is defined mod 2xpi
  93 // If two coordinates are separated by more than pi in phi one has to be shifted by +/- 2 pi
  94
  95             Double_t dx = mx[i] - x[j];
  96             if (dx >  TMath::Pi()) xx += 2. * TMath::Pi();
  97             if (dx < -TMath::Pi()) xx -= 2. * TMath::Pi();
  98             mx[i] = mx[i] * rk[i] + r[j][i] * xx;
  99             my[i] = my[i] * rk[i] + r[j][i] * y[j];
 100             rk[i] += r[j][i];
 101             mx[i] /= rk[i];
 102             my[i] /= rk[i];
 103             if (mx[i] > 2. * TMath::Pi()) mx[i] -= 2. * TMath::Pi();
 104             if (mx[i] < 0.              ) mx[i] += 2. * TMath::Pi();
 105         } // Data
 106         di += d(mx[i], my[i], oldx, oldy);
 107       } // means
 108         //
 109         // ending condition
 110       if (di < 1.e-8 || nit > 1000) break;
 111     } // while
 112
 113 // Clean-up
 114     delete[] nr;
 115     for (j = 0; j < n; j++) delete[] r[j];
 116     delete[] r;
 117 //
 118     return (nit < 1000);
 119
 120 }
 121
 122 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 )
 123 {
 124     //
 125     // The soft K-means algorithm
 126     //
 127     Int_t i,j;
 128     //
 129     // (1) Initialisation of the k means
 130     //
 131     // (there is an optimized version for initialisation called kmeans++)
 132     for (i = 0; i < k; i++) {
 133         mx[i] = 2. * TMath::Pi() * gRandom->Rndm();
 134         my[i] = -1. + 2. * gRandom->Rndm();
 135     }
 136
 137     //
 138     // (2a) The responsibilities
 139     Double_t** r = new Double_t*[n]; // responsibilities
 140     for (j = 0; j < n; j++) {r[j] = new Double_t[k];}
 141     //
 142     // (2b) Normalisation
 143     Double_t* nr = new Double_t[n];
 144     //
 145     // (2c) Weights
 146     Double_t* pi = new Double_t[k];
 147     //
 148     // (2d) Initialise the responsibilties and weights
 149     for (j = 0; j < n; j++) {
 150       nr[j] = 0.;
 151       for (i = 0; i < k; i++) {
 152         r[j][i] = TMath::Exp(- fBeta * d(mx[i], my[i], x[j], y[j]));
 153         nr[j] += r[j][i];
 154       } // mean i
 155     } // data point j
 156
 157     for (i = 0; i < k; i++) {
 158       rk[i]    = 0.;
 159       sigma2[i] = 1./fBeta;
 160
 161       for (j = 0; j < n; j++) {
 162         r[j][i] /=  nr[j];
 163         rk[i] += r[j][i];
 164       } // mean i
 165       pi[i] = rk[i] / Double_t(n);
 166     } // data point j
 167     // (3) Iterations
 168     Int_t nit = 0;
 169
 170     while(1) {
 171         nit++;
 172       //
 173       // Assignment step
 174       //
 175       for (j = 0; j < n; j++) {
 176         nr[j] = 0.;
 177         for (i = 0; i < k; i++) {
 178           r[j][i] = pi[i] * TMath::Exp(- d(mx[i], my[i], x[j], y[j]) / sigma2[i] )
 179             / (2. * sigma2[i] * TMath::Pi() * TMath::Pi());
 180           nr[j] += r[j][i];
 181         } // mean i
 182       } // data point j
 183
 184       for (i = 0; i < k; i++) {
 185         for (j = 0; j < n; j++) {
 186           r[j][i] /=  nr[j];
 187         } // mean i
 188       } // data point j
 189
 190         //
 191         // Update step
 192       Double_t di = 0;
 193
 194       for (i = 0; i < k; i++) {
 195           Double_t oldx = mx[i];
 196           Double_t oldy = my[i];
 197
 198           mx[i] = x[0];
 199           my[i] = y[0];
 200           rk[i] = r[0][i];
 201
 202         for (j = 1; j < n; j++) {
 203             Double_t xx =  x[j];
 204 //
 205 // Here we have to take into acount the cylinder topology where phi is defined mod 2xpi
 206 // If two coordinates are separated by more than pi in phi one has to be shifted by +/- 2 pi
 207
 208             Double_t dx = mx[i] - x[j];
 209             if (dx >  TMath::Pi()) xx += 2. * TMath::Pi();
 210             if (dx < -TMath::Pi()) xx -= 2. * TMath::Pi();
 211             mx[i] = mx[i] * rk[i] + r[j][i] * xx;
 212             my[i] = my[i] * rk[i] + r[j][i] * y[j];
 213             rk[i] += r[j][i];
 214             mx[i] /= rk[i];
 215             my[i] /= rk[i];
 216
 217             if (mx[i] > 2. * TMath::Pi()) mx[i] -= 2. * TMath::Pi();
 218             if (mx[i] < 0.              ) mx[i] += 2. * TMath::Pi();
 219         } // Data
 220         di += d(mx[i], my[i], oldx, oldy);
 221
 222       } // means
 223       //
 224       // Sigma
 225       for (i = 0; i < k; i++) {
 226         sigma2[i] = 0.;
 227         for (j = 1; j < n; j++) {
 228           sigma2[i] += r[j][i] * d(mx[i], my[i], x[j], y[j]);
 229         } // Data
 230         sigma2[i] /= rk[i];
 231         if (sigma2[i] < 0.0025) sigma2[i] = 0.0025;
 232       } // Clusters
 233       //
 234       // Fractions
 235       for (i = 0; i < k; i++) pi[i] = rk[i] / Double_t(n);
 236       //
 237 // ending condition
 238       if (di < 1.e-8 || nit > 1000) break;
 239     } // while
 240
 241 // Clean-up
 242     delete[] nr;
 243     delete[] pi;
 244     for (j = 0; j < n; j++) delete[] r[j];
 245     delete[] r;
 246 //
 247     return (nit < 1000);
 248 }
 249
 250 Double_t AliKMeansClustering::d(Double_t mx, Double_t my, Double_t x, Double_t y)
 251 {
 252     //
 253     // Distance definition
 254     // Quasi - Euclidian on the eta-phi cylinder
 255
 256     Double_t dx = TMath::Abs(mx-x);
 257     if (dx > TMath::Pi()) dx = 2. * TMath::Pi() - dx;
 258
 259     return (0.5*(dx * dx + (my - y) * (my - y)));
 260 }