5 * Created by Pat Schloss on 8/14/08.
6 * Copyright 2008 Patrick D. Schloss. All rights reserved.
10 #include "cluster.hpp"
11 #include "rabundvector.hpp"
12 #include "listvector.hpp"
13 #include "sparsematrix.hpp"
15 /***********************************************************************/
17 Cluster::Cluster(RAbundVector* rav, ListVector* lv, SparseMatrix* dm, float c, string f) :
18 rabund(rav), list(lv), dMatrix(dm), method(f)
21 cout << "sizeof(MatData): " << sizeof(MatData) << endl;
22 cout << "sizeof(PCell*): " << sizeof(PCell*) << endl;
24 int nCells = dMatrix->getNNodes();
25 time_t start = time(NULL);
27 MatVec matvec = MatVec(nCells);
29 for (MatData currentCell = dMatrix->begin(); currentCell != dMatrix->end(); currentCell++) {
30 matvec[i++] = currentCell;
32 for (i= matvec.size();i>0;i--) {
33 dMatrix->rmCell(matvec[i-1]);
35 MatData it = dMatrix->begin();
36 while (it != dMatrix->end()) {
37 it = dMatrix->rmCell(it);
39 cout << "Time to remove " << nCells << " cells: " << time(NULL) - start << " seconds" << endl;
41 MatData it = dMatrix->begin();
42 cout << it->row << "/" << it->column << "/" << it->dist << endl;
43 dMatrix->rmCell(dMatrix->begin());
44 cout << it->row << "/" << it->column << "/" << it->dist << endl;
48 // Create a data structure to quickly access the PCell information
49 // for a certain sequence. It consists of a vector of lists, where
50 // a list contains pointers (iterators) to the all distances related
51 // to a certain sequence. The Vector is accessed via the index of a
52 // sequence in the distance matrix.
54 //string temp = "temp";
55 //openOutputFile(temp, outtemp);
56 //cout << lv->size() << endl;
57 seqVec = vector<MatVec>(lv->size());
58 for (MatData currentCell = dMatrix->begin(); currentCell != dMatrix->end(); currentCell++) {
59 //outtemp << currentCell->row << '\t' << currentCell->column << '\t' << currentCell->dist << endl;
60 seqVec[currentCell->row].push_back(currentCell);
61 seqVec[currentCell->column].push_back(currentCell);
64 mapWanted = false; //set to true by mgcluster to speed up overlap merge
66 //save so you can modify as it changes in average neighbor
68 m = MothurOut::getInstance();
71 /***********************************************************************/
73 void Cluster::getRowColCells() {
75 PCell* smallCell = dMatrix->getSmallestCell(); //find the smallest cell - this routine should probably not be in the SpMat class
77 smallRow = smallCell->row; // get its row
78 smallCol = smallCell->column; // get its column
79 smallDist = smallCell->dist; // get the smallest distance
80 //cout << "small row = " << smallRow << "small col = " << smallCol << "small dist = " << smallDist << endl;
82 rowCells = seqVec[smallRow]; // all distances related to the row index
83 colCells = seqVec[smallCol]; // all distances related to the column index
84 nRowCells = rowCells.size();
85 nColCells = colCells.size();
86 //cout << "num rows = " << nRowCells << "num col = " << nColCells << endl;
88 //for (int i = 0; i < nColCells; i++) { cout << colCells[i]->row << '\t' << colCells[i]->column << endl; }
89 //for (int i = 0; i < nRowCells; i++) { cout << rowCells[i]->row << '\t' << rowCells[i]->column << endl; }
92 m->errorOut(e, "Cluster", "getRowColCells");
97 /***********************************************************************/
98 // Remove the specified cell from the seqVec and from the sparse
100 void Cluster::removeCell(const MatData& cell, int vrow, int vcol, bool rmMatrix){
103 ull drow = cell->row;
104 ull dcol = cell->column;
105 if (((vrow >=0) && (drow != smallRow)) ||
106 ((vcol >=0) && (dcol != smallCol))) {
116 nCells = seqVec[drow].size();
117 for (vrow=0; vrow<nCells;vrow++) {
118 crow = seqVec[drow][vrow]->row;
119 ccol = seqVec[drow][vrow]->column;
120 if (((crow == drow) && (ccol == dcol)) ||
121 ((ccol == drow) && (crow == dcol))) {
127 seqVec[drow].erase(seqVec[drow].begin()+vrow);
129 nCells = seqVec[dcol].size();
130 for (vcol=0; vcol<nCells;vcol++) {
131 crow = seqVec[dcol][vcol]->row;
132 ccol = seqVec[dcol][vcol]->column;
133 if (((crow == drow) && (ccol == dcol)) ||
134 ((ccol == drow) && (crow == dcol))) {
140 seqVec[dcol].erase(seqVec[dcol].begin()+vcol);
143 //cout << " removing = " << cell->row << '\t' << cell->column << '\t' << cell->dist << endl;
144 dMatrix->rmCell(cell);
145 // cout << "done" << endl;
149 catch(exception& e) {
150 m->errorOut(e, "Cluster", "removeCell");
154 /***********************************************************************/
156 void Cluster::clusterBins(){
158 // cout << smallCol << '\t' << smallRow << '\t' << smallDist << '\t' << rabund->get(smallRow) << '\t' << rabund->get(smallCol);
160 rabund->set(smallCol, rabund->get(smallRow)+rabund->get(smallCol));
161 rabund->set(smallRow, 0);
162 rabund->setLabel(toString(smallDist));
164 // cout << '\t' << rabund->get(smallRow) << '\t' << rabund->get(smallCol) << endl;
166 catch(exception& e) {
167 m->errorOut(e, "Cluster", "clusterBins");
174 /***********************************************************************/
176 void Cluster::clusterNames(){
178 // cout << smallCol << '\t' << smallRow << '\t' << smallDist << '\t' << list->get(smallRow) << '\t' << list->get(smallCol);
179 if (mapWanted) { updateMap(); }
181 list->set(smallCol, list->get(smallRow)+','+list->get(smallCol));
182 list->set(smallRow, "");
183 list->setLabel(toString(smallDist));
185 // cout << '\t' << list->get(smallRow) << '\t' << list->get(smallCol) << endl;
187 catch(exception& e) {
188 m->errorOut(e, "Cluster", "clusterNames");
194 /***********************************************************************/
195 //This function clusters based on the method of the derived class
196 //At the moment only average and complete linkage are covered, because
197 //single linkage uses a different approach.
198 void Cluster::update(double& cutOFF){
201 //cout << "got rowcells" << endl;
203 vector<int> foundCol(nColCells, 0);
208 // The vector has to be traversed in reverse order to preserve the index
209 // for faster removal in removeCell()
210 for (int i=nRowCells-1;i>=0;i--) {
211 if (!((rowCells[i]->row == smallRow) && (rowCells[i]->column == smallCol))) {
212 if (rowCells[i]->row == smallRow) {
213 search = rowCells[i]->column;
215 search = rowCells[i]->row;
219 for (int j=0;j<nColCells;j++) {
220 if (!((colCells[j]->row == smallRow) && (colCells[j]->column == smallCol))) { //if you are not hte smallest distance
221 if (colCells[j]->row == search || colCells[j]->column == search) {
224 changed = updateDistance(colCells[j], rowCells[i]);
225 // If the cell's distance changed and it had the same distance as
226 // the smallest distance, invalidate the mins vector in SparseMatrix
228 if (colCells[j]->vectorMap != NULL) {
229 *(colCells[j]->vectorMap) = NULL;
230 colCells[j]->vectorMap = NULL;
237 //if not merged it you need it for warning
238 if ((!merged) && (method == "average")) {
239 //m->mothurOut("Warning: trying to merge cell " + toString(rowCells[i]->row+1) + " " + toString(rowCells[i]->column+1) + " distance " + toString(rowCells[i]->dist) + " with value above cutoff. Results may vary from using cutoff at cluster command instead of read.dist."); m->mothurOutEndLine();
240 if (cutOFF > rowCells[i]->dist) {
241 cutOFF = rowCells[i]->dist;
242 //m->mothurOut("changing cutoff to " + toString(cutOFF)); m->mothurOutEndLine();
246 removeCell(rowCells[i], i , -1);
253 // Special handling for singlelinkage case, not sure whether this
255 for (int i=nColCells-1;i>=0;i--) {
256 if (foundCol[i] == 0) {
257 if (method == "average") {
258 if (!((colCells[i]->row == smallRow) && (colCells[i]->column == smallCol))) {
259 //m->mothurOut("Warning: merging cell " + toString(colCells[i]->row+1) + " " + toString(colCells[i]->column+1) + " distance " + toString(colCells[i]->dist) + " value above cutoff. Results may vary from using cutoff at cluster command instead of read.dist."); m->mothurOutEndLine();
260 if (cutOFF > colCells[i]->dist) {
261 cutOFF = colCells[i]->dist;
262 //m->mothurOut("changing cutoff to " + toString(cutOFF)); m->mothurOutEndLine();
266 removeCell(colCells[i], -1, i);
270 catch(exception& e) {
271 m->errorOut(e, "Cluster", "update");
275 /***********************************************************************/
276 void Cluster::setMapWanted(bool f) {
281 for (int i = 0; i < list->getNumBins(); i++) {
284 string names = list->get(i);
285 while (names.find_first_of(',') != -1) {
287 string name = names.substr(0,names.find_first_of(','));
288 //save name and bin number
290 names = names.substr(names.find_first_of(',')+1, names.length());
298 catch(exception& e) {
299 m->errorOut(e, "Cluster", "setMapWanted");
303 /***********************************************************************/
304 void Cluster::updateMap() {
306 //update location of seqs in smallRow since they move to smallCol now
307 string names = list->get(smallRow);
308 while (names.find_first_of(',') != -1) {
310 string name = names.substr(0,names.find_first_of(','));
311 //save name and bin number
312 seq2Bin[name] = smallCol;
313 names = names.substr(names.find_first_of(',')+1, names.length());
317 seq2Bin[names] = smallCol;
320 catch(exception& e) {
321 m->errorOut(e, "Cluster", "updateMap");
325 /***********************************************************************/