+
+
+
+void AliKMeansClustering::OptimalInit(Int_t k, Int_t n, const Double_t* x, const Double_t* y, Double_t* mx, Double_t* my)
+{
+ //
+ // Optimal initialisation using the k-means++ algorithm
+ // http://en.wikipedia.org/wiki/K-means%2B%2B
+ //
+ // k-means++ is an algorithm for choosing the initial values for k-means clustering in statistics and machine learning.
+ // It was proposed in 2007 by David Arthur and Sergei Vassilvitskii as an approximation algorithm for the NP-hard k-means problem---
+ // a way of avoiding the sometimes poor clusterings found by the standard k-means algorithm.
+ //
+ //
+ TH1F d2("d2", "", n, -0.5, Float_t(n)-0.5);
+ d2.Reset();
+
+ // (1) Chose first center as a random point among the input data.
+ Int_t ir = Int_t(Float_t(n) * gRandom->Rndm());
+ mx[0] = x[ir];
+ my[0] = y[ir];
+
+ // (2) Iterate
+ Int_t icl = 1;
+ while(icl < k)
+ {
+ // find min distance to existing clusters
+ for (Int_t j = 0; j < n; j++) {
+ Double_t dmin = 1.e10;
+ for (Int_t i = 0; i < icl; i++) {
+ Double_t dij = d(mx[i], my[i], x[j], y[j]);
+ if (dij < dmin) dmin = dij;
+ } // clusters
+ d2.Fill(Float_t(j), dmin);
+ } // data points
+ // select a new cluster from data points with probability ~d2
+ ir = Int_t(d2.GetRandom() + 0.5);
+ mx[icl] = x[ir];
+ my[icl] = y[ir];
+ icl++;
+ } // icl
+}
+
+
+ClassImp(AliKMeansResult)
+
+
+
+AliKMeansResult::AliKMeansResult(Int_t k):
+ TObject(),
+ fK(k),
+ fMx (new Double_t[k]),
+ fMy (new Double_t[k]),
+ fSigma2(new Double_t[k]),
+ fRk (new Double_t[k]),
+ fTarget(new Double_t[k]),
+ fInd (new Int_t[k])
+{
+// Constructor
+}
+
+AliKMeansResult::AliKMeansResult(const AliKMeansResult &res):
+ TObject(res),
+ fK(res.GetK()),
+ fMx(new Double_t[res.GetK()]),
+ fMy(new Double_t[res.GetK()]),
+ fSigma2(new Double_t[res.GetK()]),
+ fRk(new Double_t[res.GetK()]),
+ fTarget(new Double_t[res.GetK()]),
+ fInd(new Int_t[res.GetK()])
+{
+ // Copy constructor
+ for (Int_t i = 0; i <fK; i++) {
+ fMx[i] = (res.GetMx()) [i];
+ fMy[i] = (res.GetMy()) [i];
+ fSigma2[i] = (res.GetSigma2())[i];
+ fRk[i] = (res.GetRk()) [i];
+ fTarget[i] = (res.GetTarget())[i];
+ fInd[i] = (res.GetInd()) [i];
+ }
+}
+
+AliKMeansResult& AliKMeansResult::operator=(const AliKMeansResult& res)
+{
+ //
+ // Assignment operator
+ if (this != &res) {
+ TObject::operator=(res);
+ if (fK != res.fK) {
+ delete [] fMx;
+ delete [] fMy;
+ delete [] fSigma2;
+ delete [] fRk;
+ delete [] fTarget;
+ delete [] fInd;
+ fK = res.fK;
+ fMx = new Double_t[fK];
+ fMy = new Double_t[fK];
+ fSigma2 = new Double_t[fK];
+ fRk = new Double_t[fK];
+ fTarget = new Double_t[fK];
+ fInd = new Int_t[fK];
+ }
+
+ fK = res.fK;
+ memcpy(fMx, res.fMx, fK*sizeof(Double_t));
+ memcpy(fMy, res.fMy, fK*sizeof(Double_t));
+ memcpy(fSigma2, res.fSigma2, fK*sizeof(Double_t));
+ memcpy(fRk, res.fRk, fK*sizeof(Double_t));
+ memcpy(fTarget, res.fTarget, fK*sizeof(Double_t));
+ memcpy(fInd, res.fInd, fK*sizeof(Int_t));
+ }
+ return *this;
+}
+
+
+AliKMeansResult::~AliKMeansResult()
+{
+// Destructor
+ delete[] fMx;
+ delete[] fMy;
+ delete[] fSigma2;
+ delete[] fRk;
+ delete[] fInd;
+ delete[] fTarget;
+}
+
+void AliKMeansResult::Sort()
+{
+ // Build target array and sort
+ // Sort clusters
+ for (Int_t i = 0; i < fK; i++) {
+ if (fRk[i] > 2.9) {
+ fTarget[i] = fRk[i] / fSigma2[i];
+ }
+ else fTarget[i] = 0.;
+ }
+
+ TMath::Sort(fK, fTarget, fInd);
+}
+
+void AliKMeansResult::Sort(Int_t n, const Double_t* x, const Double_t* y)
+{
+ // Build target array and sort
+ for (Int_t i = 0; i < fK; i++)
+ {
+ Int_t nc = 0;
+ for (Int_t j = 0; j < n; j++)
+ {
+ if (2. * AliKMeansClustering::d(fMx[i], fMy[i], x[j], y[j]) < 2.28 * fSigma2[i]) nc++;
+ }
+
+ if (nc > 2) {
+ fTarget[i] = Double_t(nc) / (2.28 * fSigma2[i]);
+ } else {
+ fTarget[i] = 0.;
+ }
+ }
+
+ TMath::Sort(fK, fTarget, fInd);
+}
+
+void AliKMeansResult::CopyResults(const AliKMeansResult* res)
+{
+ fK = res->GetK();
+ for (Int_t i = 0; i <fK; i++) {
+ fMx[i] = (res->GetMx()) [i];
+ fMy[i] = (res->GetMy()) [i];
+ fSigma2[i] = (res->GetSigma2())[i];
+ fRk[i] = (res->GetRk()) [i];
+ fTarget[i] = (res->GetTarget())[i];
+ fInd[i] = (res->GetInd()) [i];
+ }
+}