-vector<double> Pam::calcSilhouettes(vector<vector<double> > dists) {
- try {
- vector<double> silhouettes; silhouettes.resize(numSamples, 0.0);
- if (numPartitions < 2) { return silhouettes; }
-
- vector<int> whoDp;
- for (int i = 0; i < numSamples; i++) { // assumes at least 2 partitions
- vector<seqDist> tempMeds;
- for (set<int>::iterator it = medoids.begin(); it != medoids.end(); it++) {
- if (m->control_pressed) { break; }
- seqDist temp(*it, *it, dists[i][*it]); //medoid, medoid, distance from medoid to sample
- tempMeds.push_back(temp);
- }
- sort(tempMeds.begin(), tempMeds.end(), compareSequenceDistance); //sort by distance
- whoDp.push_back(tempMeds[1].seq1); //second closest medoid
- }
-
-
- vector<double> a, b; a.resize(numSamples, 0.0); b.resize(numSamples, 0.0);
-
- //calc a - all a[i] are the same in the same partition
- for (int k = 0; k < numPartitions; k++) {
- if (m->control_pressed) { break; }
-
- int count = 0;
- double totalA = 0.0;
- for (int i = 0; i < numSamples; i++) {
- for (int j = 0; j < numSamples; j++) {
- if (m->control_pressed) { break; }
- if ((zMatrix[k][i] != 0) && (zMatrix[k][j] != 0)){ //are both samples in the partition, if so add there distance
- totalA += dists[j][i]; //distance from this sample to the other samples in the partition
- count++;
- }
- }
- }
- totalA /= (double) count;
-
- //set a[i] to average for cluster
- for (int i = 0; i < numSamples; i++) {
- if (zMatrix[k][i] != 0) { a[i] = totalA; }
- }
- }
-
- //calc b
- for (int i = 0; i < numSamples; i++) {
- if (m->control_pressed) { break; }
-
- int nextClosestMedoid = whoDp[i];
- int thisPartition = medoid2Partition[nextClosestMedoid];
- int count = 0;
- double totalB = 0.0;
- for (int j = 0; j < numSamples; j++) {
- if (zMatrix[thisPartition][j] != 0) { //this sample is in this partition
- totalB += dists[i][j];
- count++;
- }
- }
- b[i] = totalB / (double) count;
- }
-
- //calc silhouettes
- for (int i = 0; i < numSamples; i++) {
- if (m->control_pressed) { break; }
-
- double denom = a[i];
- if (b[i] > denom) { denom = b[i]; } //max(a[i],b[i])
-
- silhouettes[i] = (b[i] - a[i]) / denom;
-
- //cout << "silhouettes " << i << '\t' << silhouettes[i] << endl;
- }
-
- return silhouettes;
- }
- catch(exception& e) {
- m->errorOut(e, "Pam", "calcSilhouettes");
- exit(1);
- }
-}