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 //build and swap functions based on pam.c by maechler from R cluster package
  39 //sets Dp[0] does not set Dp[1]. chooses intial medoids.
  40 int Pam::buildPhase() {
  41     try {
  42
  43         if (m->debug) { m->mothurOut("[DEBUG]: building medoids\n"); }
  44
  45         vector<double> gains; gains.resize(numSamples);
  46
  47         largestDist *= 1.1 + 1; //make this distance larger than any distance in the matrix
  48         Dp.resize(numSamples);
  49         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
  50
  51         zMatrix.resize(numPartitions);
  52         for(int i=0;i<numPartitions;i++){
  53             zMatrix[i].assign(numSamples, 0);
  54         }
  55
  56         for (int k = 0; k < numPartitions; k++) {
  57
  58             int medoid = -1;
  59             double totalGain = 0.0;
  60             double clusterGain = 0.0;
  61
  62             for (int i = 0; i < numSamples; i++) {  //does this need to be square?? can we do lt?
  63                 if (m->control_pressed) { break; }
  64
  65                 if (medoids.count(i) == 0) { //is this sample is NOT a medoid?
  66                     gains[i] = 0.0;
  67
  68                     for (int j = 0; j < numSamples; j++) {
  69                         totalGain = Dp[j][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 += (- Dp[j][0]+ small);
 132                                 }else if (dists[h][j] < Dp[j][0]) {
 133                                     dz += (- Dp[j][0] + dists[h][j]);
 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
 226         return 0;
 227     }
 228     catch(exception& e) {
 229         m->errorOut(e, "Pam", "updateDp");
 230         exit(1);
 231     }
 232 }
 233
 234 /**************************************************************************************************/
 235 /*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:
 236
 237  CHk=Bk/(k−1)/Wk/(n−k)
 238
 239  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.*/
 240 //based on R index.G1.r function
 241 double Pam::calcCHIndex(vector< vector<double> > dists){ //countMatrix = [numSamples][numOtus]
 242     try {
 243         double CH = 0.0;
 244
 245         if (numPartitions < 2) { return CH; }
 246
 247         map<int, int> clusterMap; //map sample to partition
 248         for (int i = 0; i < numPartitions; i++) {
 249             for (int j = 0; j < numSamples; j++) {
 250                 if (m->control_pressed) { return 0.0; }
 251                 if (zMatrix[i][j] != 0) { clusterMap[j] = i; }
 252             }
 253         }
 254
 255         //make countMatrix a relabund
 256         vector<vector<double> > relativeAbundance(numSamples); //[numSamples][numOTUs]
 257         //get relative abundance
 258         for(int i=0;i<numSamples;i++){
 259             if (m->control_pressed) {  return 0; }
 260             int groupTotal = 0;
 261
 262             relativeAbundance[i].assign(numOTUs, 0.0);
 263
 264             for(int j=0;j<numOTUs;j++){
 265                 groupTotal += countMatrix[i][j];
 266             }
 267             for(int j=0;j<numOTUs;j++){
 268                 relativeAbundance[i][j] = countMatrix[i][j] / (double)groupTotal;
 269             }
 270         }
 271
 272         //find centers
 273         vector<vector<double> > centers; centers.resize(numPartitions);
 274         int countPartitions = 0;
 275         for (set<int>::iterator it = medoids.begin(); it != medoids.end(); it++) {
 276             for (int j = 0; j < numOTUs; j++) {
 277                 centers[countPartitions].push_back(relativeAbundance[*it][j]); //save the relative abundance of the medoid for this partition for this OTU
 278             }
 279             countPartitions++;
 280         }
 281
 282         //centers.clear();
 283         //centers = calcCenters(dists, clusterMap, relativeAbundance);
 284
 285         double allMeanDist = rMedoid(relativeAbundance, dists);
 286
 287         if (m->debug) { m->mothurOut("[DEBUG]: allMeandDist = " + toString(allMeanDist) + "\n"); }
 288
 289         for (int i = 0; i < relativeAbundance.size(); i++) {//numSamples
 290             for (int j = 0; j < relativeAbundance[i].size(); j++) { //numOtus
 291                 if (m->control_pressed) {  return 0; }
 292                 //x <- (x - centers[cl, ])^2
 293                 relativeAbundance[i][j] = ((relativeAbundance[i][j] - centers[clusterMap[i]][j])*(relativeAbundance[i][j] - centers[clusterMap[i]][j]));
 294             }
 295         }
 296
 297         double wgss = 0.0;
 298         for (int j = 0; j < numOTUs; j++) {
 299             for(int i=0;i<numSamples;i++){
 300                 if (m->control_pressed) { return 0.0; }
 301                 wgss += relativeAbundance[i][j];
 302             }
 303         }
 304
 305         double bgss = allMeanDist - wgss;
 306
 307         CH = (bgss / (double)(numPartitions - 1)) / (wgss / (double) (numSamples - numPartitions));
 308
 309         return CH;
 310     }
 311     catch(exception& e){
 312         m->errorOut(e, "Pam", "calcCHIndex");
 313         exit(1);
 314     }
 315 }
 316
 317 /**************************************************************************************************/
 318
 319