pam.cpp

   1 //
   2 //  pam.cpp
   3 //  Mothur
   4 //
   5 //  Created by SarahsWork on 12/10/13.
   6 //  Copyright (c) 2013 Schloss Lab. All rights reserved.
   7 //
   8
   9 #include "pam.h"
  10 #define DBL_EPSILON 1e-9
  11
  12 /**************************************************************************************************/
  13 Pam::Pam(vector<vector<int> > c, vector<vector<double> > d, int p) : CommunityTypeFinder() {
  14     try {
  15         countMatrix = c;
  16         numSamples = (int)d.size();
  17         numOTUs = (int)c[0].size();
  18         numPartitions = p;
  19         dists = d;
  20
  21         largestDist = 0;
  22         for (int i = 0; i < dists.size(); i++) {
  23             for (int j = i; j < dists.size(); j++) {
  24                 if (m->control_pressed) { break; }
  25                 if (dists[i][j] > largestDist) { largestDist = dists[i][j]; }
  26             }
  27         }
  28
  29         buildPhase(); //choosing the medoids
  30         swapPhase(); //optimize clusters
  31     }
  32         catch(exception& e) {
  33                 m->errorOut(e, "Pam", "Pam");
  34                 exit(1);
  35         }
  36 }
  37 /**************************************************************************************************/
  38 //sets Dp[0] does not set Dp[1]. chooses intial medoids.
  39 int Pam::buildPhase() {
  40     try {
  41
  42         if (m->debug) { m->mothurOut("[DEBUG]: building medoids\n"); }
  43
  44         vector<double> gains; gains.resize(numSamples);
  45
  46         largestDist *= 1.1 + 1; //make this distance larger than any distance in the matrix
  47         Dp.resize(numSamples);
  48         for (int i = 0; i < numSamples; i++) { Dp[i].push_back(largestDist); Dp[i].push_back(largestDist); } //2 smallest dists for this sample in this partition
  49
  50         zMatrix.resize(numPartitions);
  51         for(int i=0;i<numPartitions;i++){
  52             zMatrix[i].assign(numSamples, 0);
  53         }
  54
  55         for (int k = 0; k < numPartitions; k++) {
  56
  57             int medoid = -1;
  58             double totalGain = 0.0;
  59             double clusterGain = 0.0;
  60
  61             for (int i = 0; i < numSamples; i++) {  //does this need to be square?? can we do lt?
  62                 if (m->control_pressed) { break; }
  63
  64                 if (medoids.count(i) == 0) { //is this sample is NOT a medoid?
  65                     gains[i] = 0.0;
  66
  67                     for (int j = 0; j < numSamples; j++) {
  68                         //cout << i << '\t' << j << '\t' <<   Dp[i][0] << '\t' << dists[i][j] << '\t' << totalGain << endl;
  69                         totalGain = Dp[i][0] - dists[i][j];
  70                         if (totalGain > 0.0) { gains[i] += totalGain; }
  71                     }
  72                     if (m->debug) { m->mothurOut("[DEBUG]: " + toString(i) +  " totalGain = " + toString(totalGain) + "\n"); }
  73
  74                     if (clusterGain <= gains[i]) {
  75                         clusterGain = gains[i];
  76                         medoid = i;
  77                     }
  78                 }
  79             }
  80
  81             //save medoid value
  82             medoids.insert(medoid);
  83
  84             if (m->debug) { m->mothurOut("[DEBUG]: new medoid " + toString(medoid) + "\n"); }
  85
  86             //update dp values
  87             for (int i = 0; i < numSamples; i++) {
  88                 if (Dp[i][0] > dists[i][medoid]) { Dp[i][0] = dists[i][medoid]; }
  89             }
  90         }
  91         if (m->debug) { m->mothurOut("[DEBUG]: done building medoids\n"); }
  92         return 0;
  93     }
  94         catch(exception& e) {
  95                 m->errorOut(e, "Pam", "buildPhase");
  96                 exit(1);
  97         }
  98 }
  99 /**************************************************************************************************/
 100 //goal to swap medoids with non-medoids to see if we can reduce the overall cost
 101 int Pam::swapPhase() {
 102     try {
 103         if (m->debug) { m->mothurOut("[DEBUG]: swapping  medoids\n"); }
 104         //calculate cost of initial choice - average distance of samples to their closest medoid
 105         double sky = 0.0;
 106         double dzsky = 1.0;
 107         for (int i = 0; i < numSamples; i++) { sky += Dp[i][0]; }  sky /= (double) numSamples;
 108
 109         bool done = false;
 110         int hbest, nbest; hbest = -1; nbest = -1;
 111         while (!done) {
 112             if (m->control_pressed) { break; }
 113
 114             updateDp();
 115
 116             dzsky = 1;
 117
 118             for (int h = 0; h < numSamples; h++) {
 119                 if (m->control_pressed) { break; }
 120                 if (medoids.count(h) == 0) { //this is NOT a medoid
 121                     for (int i = 0; i < numSamples; i++) {
 122                         if (medoids.count(i) != 0) { //this is a medoid
 123
 124                             double dz = 0.0; //Tih sum of distances between objects and closest medoid caused by swapping i and h. Basically the change in cost. If this < 0 its a "good" swap. When all Tih are > 0, then we stop the algo, because we have the optimal medoids.
 125                             for (int j = 0; j < numSamples; j++) {
 126                                 if (m->control_pressed) { break; }
 127                                 if (dists[i][j] == Dp[j][0]) {
 128                                     double small = 0.0;
 129                                     if (Dp[j][1] > dists[h][j]) {   small = dists[h][j];    }
 130                                     else                        {   small = Dp[j][1];       }
 131                                     dz += (small - Dp[j][0]);
 132                                 }else if (dists[h][j] < Dp[j][0]) {
 133                                     dz += (dists[h][j] - Dp[j][0]);
 134                                 }
 135                             }
 136                             if (dzsky > dz) {
 137                                 dzsky = dz;
 138                                 hbest = h;
 139                                 nbest = i;
 140                             }
 141                         }//end if medoid
 142                     }//end for i
 143                 }//end if NOT medoid
 144             }//end if h
 145
 146             if (dzsky < -16 *DBL_EPSILON * fabs(sky)) {
 147                 medoids.insert(hbest);
 148                 medoids.erase(nbest);
 149                 if (m->debug) { m->mothurOut("[DEBUG]: swapping " + toString(hbest) + " " + toString(nbest) + "\n"); }
 150                 sky += dzsky;
 151             }else { done = true; } //stop algo.
 152         }
 153
 154
 155         //fill zmatrix
 156         int count = 0;
 157         vector<int> tempMedoids;
 158         for (set<int>::iterator it = medoids.begin(); it != medoids.end(); it++) {
 159             medoid2Partition[*it] = count;
 160             zMatrix[count][*it] = 1; count++; //set medoid in this partition.
 161             tempMedoids.push_back(*it);
 162         }
 163
 164         //which partition do you belong to?
 165         laplace = 0;
 166         for (int i = 0; i < numSamples; i++) {
 167             int partition = 0;
 168             double dist = dists[i][tempMedoids[0]]; //assign to first medoid
 169             for (int j = 1; j < tempMedoids.size(); j++) {
 170                 if (dists[i][tempMedoids[j]] < dist) { //is this medoid closer?
 171                     dist = dists[i][tempMedoids[j]];
 172                     partition = j;
 173                 }
 174             }
 175             zMatrix[partition][i] = 1;
 176             laplace += dist;
 177         }
 178         laplace /= (double) numSamples;
 179
 180         if (m->debug) {
 181             for(int i=0;i<numPartitions;i++){
 182                 m->mothurOut("[DEBUG]: partition 1: ");
 183                 for (int j = 0; j < numSamples; j++) {
 184                      m->mothurOut(toString(zMatrix[i][j]) + " ");
 185                 }
 186                 m->mothurOut("\n");
 187             }
 188             m->mothurOut("[DEBUG]: medoids : ");
 189             for (set<int>::iterator it = medoids.begin(); it != medoids.end(); it++) { m->mothurOut(toString(*it) + " ");
 190             }
 191             m->mothurOut("\n");
 192
 193             m->mothurOut("[DEBUG]: laplace : " + toString(laplace));  m->mothurOut("\n");
 194         }
 195
 196         if (m->debug) { m->mothurOut("[DEBUG]: done swapping  medoids\n"); }
 197         return 0;
 198     }
 199     catch(exception& e) {
 200         m->errorOut(e, "Pam", "swapPhase");
 201         exit(1);
 202     }
 203 }
 204
 205 /**************************************************************************************************/
 206 int Pam::updateDp() {
 207     try {
 208         for (int j = 0; j < numSamples; j++) {
 209             if (m->control_pressed) { break; }
 210
 211             //initialize dp and ep
 212             Dp[j][0] = largestDist; Dp[j][1] = largestDist;
 213
 214             for (int i = 0; i < numSamples; i++) {
 215                 if (medoids.count(i) != 0) { //is this a medoid?
 216                     if (Dp[j][0] > dists[j][i]) {
 217                         Dp[j][0] = Dp[j][1];
 218                         Dp[j][0] = dists[j][i];
 219                     }else if (Dp[j][1] > dists[j][i]) {
 220                         Dp[j][1] = dists[j][i];
 221                     }
 222                 }
 223             }
 224         }
 225         return 0;
 226     }
 227     catch(exception& e) {
 228         m->errorOut(e, "Pam", "updateDp");
 229         exit(1);
 230     }
 231 }
 232 /**************************************************************************************************/
 233 //The silhouette width S(i)of individual data points i is calculated using the following formula:
 234 /*
 235  s(i) = b(i) - a(i)
 236  -----------
 237  max(b(i),a(i))
 238  where a(i) is the average dissimilarity (or distance) of sample i to all other samples in the same cluster, while b(i) is the average dissimilarity (or distance) to all objects in the closest other cluster.
 239
 240  The formula implies -1 =< S(i) =< 1 . A sample which is much closer to its own cluster than to any other cluster has a high S(i) value, while S(i) close to 0 implies that the given sample lies somewhere between two clusters. Large negative S(i) values indicate that the sample was assigned to the wrong cluster.
 241  */
 242
 243 vector<double> Pam::calcSilhouettes(vector<vector<double> > dists) {
 244     try {
 245         vector<double> silhouettes; silhouettes.resize(numSamples, 0.0);
 246         if (numPartitions < 2) { return silhouettes; }
 247
 248         vector<int> whoDp;
 249         for (int i = 0; i < numSamples; i++) { // assumes at least 2 partitions
 250             vector<seqDist> tempMeds;
 251             for (set<int>::iterator it = medoids.begin(); it != medoids.end(); it++) {
 252                 if (m->control_pressed) { break; }
 253                 seqDist temp(*it, *it, dists[i][*it]); //medoid, medoid, distance from medoid to sample
 254                 tempMeds.push_back(temp);
 255             }
 256             sort(tempMeds.begin(), tempMeds.end(), compareSequenceDistance); //sort by distance
 257             whoDp.push_back(tempMeds[1].seq1); //second closest medoid
 258         }
 259
 260
 261         vector<double> a, b; a.resize(numSamples, 0.0); b.resize(numSamples, 0.0);
 262
 263         //calc a - all a[i] are the same in the same partition
 264         for (int k = 0; k < numPartitions; k++) {
 265             if (m->control_pressed) { break; }
 266
 267             int count = 0;
 268             double totalA = 0.0;
 269             for (int i = 0; i < numSamples; i++) {
 270                 for (int j = 0; j < numSamples; j++) {
 271                     if (m->control_pressed) { break; }
 272                     if ((zMatrix[k][i] != 0) && (zMatrix[k][j] != 0)){ //are both samples in the partition, if so add there distance
 273                         totalA += dists[j][i]; //distance from this sample to the other samples in the partition
 274                         count++;
 275                     }
 276                 }
 277             }
 278             totalA /= (double) count;
 279
 280             //set a[i] to average for cluster
 281             for (int i = 0; i < numSamples; i++) {
 282                 if (zMatrix[k][i] != 0) { a[i] = totalA; }
 283             }
 284         }
 285
 286         //calc b
 287         for (int i = 0; i < numSamples; i++) {
 288             if (m->control_pressed) { break; }
 289
 290             int nextClosestMedoid = whoDp[i];
 291             int thisPartition = medoid2Partition[nextClosestMedoid];
 292             int count = 0;
 293             double totalB = 0.0;
 294             for (int j = 0; j < numSamples; j++) {
 295                 if (zMatrix[thisPartition][j] != 0) { //this sample is in this partition
 296                     totalB += dists[i][j];
 297                     count++;
 298                 }
 299             }
 300             b[i] = totalB / (double) count;
 301         }
 302
 303         //calc silhouettes
 304         for (int i = 0; i < numSamples; i++) {
 305             if (m->control_pressed) { break; }
 306
 307             double denom = a[i];
 308             if (b[i] > denom) { denom = b[i]; } //max(a[i],b[i])
 309
 310             silhouettes[i] = (b[i] - a[i]) / denom;
 311
 312             //cout << "silhouettes " << i << '\t' << silhouettes[i] << endl;
 313         }
 314
 315         return silhouettes;
 316     }
 317     catch(exception& e) {
 318         m->errorOut(e, "Pam", "calcSilhouettes");
 319         exit(1);
 320     }
 321 }
 322 /**************************************************************************************************/
 323 /*To assess the optimal number of clusters our dataset was most robustly partitioned into, we used the Calinski-Harabasz (CH) Index that has shown good performance in recovering the number of clusters. It is defined as:
 324
 325  CHk=Bk/(k−1)/Wk/(n−k)
 326
 327  where Bk is the between-cluster sum of squares (i.e. the squared distances between all points i and j, for which i and j are not in the same cluster) and Wk is the within-clusters sum of squares (i.e. the squared distances between all points i and j, for which i and j are in the same cluster). This measure implements the idea that the clustering is more robust when between-cluster distances are substantially larger than within-cluster distances. Consequently, we chose the number of clusters k such that CHk was maximal.*/
 328 double Pam::calcCHIndex(vector< vector<double> > dists){ //countMatrix = [numSamples][numOtus]
 329     try {
 330         double CH = 0.0;
 331
 332         if (numPartitions < 2) { return CH; }
 333
 334         map<int, int> clusterMap; //map sample to partition
 335         for (int i = 0; i < numPartitions; i++) {
 336             for (int j = 0; j < numSamples; j++) {
 337                 if (m->control_pressed) { return 0.0; }
 338                 if (zMatrix[i][j] != 0) { clusterMap[j] = i; }
 339             }
 340         }
 341
 342         double sumBetweenCluster = 0.0;
 343         double sumWithinClusters = 0.0;
 344
 345         for (int i = 0; i < numSamples; i++) { //lt
 346             for (int j = 0; j < i; j++) {
 347                 if (m->control_pressed) { return 0.0; }
 348                 int partitionI = clusterMap[i];
 349                 int partitionJ = clusterMap[j];
 350
 351                 if (partitionI == partitionJ) { //they are from the same cluster so this distance is added to Wk
 352                     sumWithinClusters += (dists[i][j] * dists[i][j]);
 353                 }else { //they are NOT from the same cluster so this distance is added to Bk
 354                     sumBetweenCluster += (dists[i][j] * dists[i][j]);
 355                 }
 356             }
 357         }
 358         //cout << numPartitions << '\t' << sumWithinClusters << '\t' << sumBetweenCluster << '\t' <<  (numSamples - numPartitions) << endl;
 359
 360         CH = (sumBetweenCluster / (double)(numPartitions - 1)) / (sumWithinClusters / (double) (numSamples - numPartitions));
 361
 362         return CH;
 363     }
 364     catch(exception& e){
 365         m->errorOut(e, "Pam", "calcCHIndex");
 366         exit(1);
 367     }
 368 }
 369
 370 /**************************************************************************************************/
 371
 372
 373