1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // Implemenatation of the K-Means Clustering Algorithm
17 // See http://en.wikipedia.org/wiki/K-means_clustering and references therein.
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.
22 // Author: Andreas Morsch (CERN)
23 // andreas.morsch@cern.ch
25 #include "AliKMeansClustering.h"
29 ClassImp(AliKMeansClustering)
31 Double_t AliKMeansClustering::fBeta = 10.;
33 void AliKMeansClustering::SoftKMeans(Int_t k, Int_t n, Double_t* x, Double_t* y, Double_t* mx, Double_t* my , Double_t* rk )
36 // The soft K-means algorithm
40 // (1) Initialisation of the k means
42 for (i = 0; i < k; i++) {
43 mx[i] = 2. * TMath::Pi() * gRandom->Rndm();
44 my[i] = -1. + 2. * gRandom->Rndm();
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];}
53 Double_t* nr = new Double_t[n];
63 for (j = 0; j < n; j++) {
65 for (i = 0; i < k; i++) {
66 r[j][i] = TMath::Exp(- fBeta * d(mx[i], my[i], x[j], y[j]));
71 for (j = 0; j < n; j++) {
72 for (i = 0; i < k; i++) {
81 for (i = 0; i < k; i++) {
82 Double_t oldx = mx[i];
83 Double_t oldy = my[i];
89 for (j = 1; j < n; j++) {
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
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];
103 if (mx[i] > 2. * TMath::Pi()) mx[i] -= 2. * TMath::Pi();
104 if (mx[i] < 0. ) mx[i] += 2. * TMath::Pi();
106 di += d(mx[i], my[i], oldx, oldy);
110 if (di < 1.e-8 || nit > 1000) break;
115 for (j = 0; j < n; j++) delete[] r[j];
119 Double_t AliKMeansClustering::d(Double_t mx, Double_t my, Double_t x, Double_t y)
122 // Distance definition
123 // Quasi - Euclidian on the eta-phi cylinder
125 Double_t dx = TMath::Abs(mx-x);
126 if (dx > TMath::Pi()) dx = 2. * TMath::Pi() - dx;
128 return (dx * dx + (my - y) * (my - y));